离散哈特莱变换(DHT)及快速哈特莱变换(FHT)学习

离散哈特莱变换(DHT)及快速哈特莱变换(FHT)学习

---

说在前边

最近复习$DSP$的时候，发现了一个号称专门针对离散实序列的变换，经分析总运算量为普通$FFT$的几乎一半，而且完全没有复数。这么强的吗？于是花了一个下午，去学习了一下。。。于是去图书馆翻了几乎所有的$dsp$课本。。。发现了这本书 西安电子科技大学出版社《数字信号处理》第二版！竟然花了一节在讲$DHT$和$FHT$！竟然还有附录FORTRAN代码（虽然一堆错）！这里把这个东西稍微普及一下。。。不过估计没人用的上

---

离散哈特莱变换(DHT)

引入

对于序列$x(n)$其$N$点$DFT$具有简单的共轭对称性，即 $$X(N - k) = X^*(k), k = 0, 1 ... ,N-1$$
所以只要计算$X(k)$的前$N/2$个值，则后$N/2$可以通过上式求得，$X(k)$的$N/2$个复数正好对应$N$个实数数据($N$点实序列$DFT$，也可以通过把$N$个实数，压成$N/2$个复数，再利用共轭对称性求解，这里不做详解)。由此可见，一个$N$点实序列的$DFT$，完全可由$N$个实数数据确定。由此，我们引出一种直接对于实序列进行实数域变换的离散哈特莱变换($DHT$)。$DHT$将$x(n)$看成实数数据，而不是像$DFT$一样看作虚部为$0$的复数，因此节省了一半的空间，运算效率也提升了近一倍，且$DFT$与$DHT$之间存在简单的关系，容易实现相互转换。

预备知识

1. $DFT$的意义
2. $FFT$实现
3. 最好学过$DSP$

定义

设$x(n), n = 0, 1,..., N-1$，为一实序列，其$DHT$定义为

$$X_H(k) = DHT[ x(n)] = \sum _{n=0}^{N-1} x(n)cas(\frac{2\pi}{N}kn), k = 0,1, ... , N-1$$

式中$cas(a) = cos(a) + sin(a)$逆变换($IDHT$)为

$$x(n) = DHT[ X_H(k)] = \frac{1}{N}\sum _{n=0}^{N-1} X(k)cas(\frac{2\pi}{N}kn), n = 0,1, ... , N-1$$

$DHT$的正交证明：

$$\sum_{k=0}^{N-1} cas(\frac{2\pi}{N}kn)cas(\frac{2\pi}{N}km) = \sum_{k=0}^{N-1}[cos(\frac{2\pi}{N}k(n-m)) + sin(\frac{2\pi}{N}k(n+m))] = \left\{\begin{array}{cc} N, & k = 0\\ 0, & k \neq 0 \\ \end{array}\right.$$

DHT与DFT的关系

用$X(k)$表示实序列$x(n)$的$DFT$，用$X_H(k)$表示$x(n)$的$DHT$，分别用$X_{He}(k)$， $X_{Ho}(k)$表示$X_H(k)$的偶对称分量与奇对称分量，即

$$X_H(k) = X_{He}(k) + X_{Ho}(k)$$

其中

$$X_{He}(k) = \frac{1}{2}[X_H(k) + X_H(N-k)]\\ X_{Ho}(k) = \frac{1}{2}[X_H(k) - X_H(N-k)]\\ X(k) = X_{He}(k) - jX_{Ho}(k)\\ X_H(k) = Re[X(k)] - Im[X(k)]\\ X(k) = \frac{1}{2}[X_H(k)+X_H(N-k)] - \frac{1}{2}j[X_H(k)-X_H(N-k)]$$

利用上面这些性质，我们可以很容易的将$DFT$与$DHT$进行变换，并且又因为这个原因$X_H(k)$的周期为$N$(隐含周期性)。

DHT的优点

1. $DHT$为实值，避免了复数运算
2. $DHT$正反变换形式基本一致
3. $DHT$与$DFT$的转换容易实现

DHT的性质

1. 线性性
2. $DHT$不改变$x(n)$奇偶性
3. 循环卷积定理

若$x_1(n) \leftrightarrow X_{1H}(k), ~~x_2(n) \leftrightarrow X_{2H}(k)$则

$$x_1(n) \otimes x_2(n) \leftrightarrow X_{2H}(k)X_{1He}(k) + X_{2H}(N-k)X_{1Ho}(k)$$

或

$$x_1(n) \otimes x_1(n) \leftrightarrow X_{1H}(k)X_{2He}(k) + X_{1H}(N-k)X_{2Ho}(k)$$

值得注意的是，相比于$DFT$的卷积定理，$DHT$的运算更加复杂，这一点把这个算法的在竞赛中的实用性大大降低了。。。毕竟我们就是为了加速卷积，不过他空间上的优势还是很名明显的。

---

快速哈特莱变换(FHT)

基2 DIT-FHT算法

与$FFT$算法一样，$FHT$也可以，用基$4$,基$8$，分裂基等方式优化，这里我们讨论最基础的基$2$ 快速$DHT$算法。
$x(n)$的$N=2^M$点$DHT$：

$$X_H(k) = \sum_{n=0}^{N-1}x(n)cas(\frac{2\pi}{N}kn)$$

对$x(n)$进行奇偶抽取

$$x_0(r) = x(2r)\\ x_1(r) = x(2r+1)$$

带入后，与$DFT$比较可得，$cas(\frac{2\pi}{N}k(2r+1))$不是指数函数，所以我们通过

$$cas(a + b) = cas(a)cas(b) + cas(-a)sin(b)$$

推导可得：

$$X_H(k) = \sum_{r=0}^{\frac{N}{2}-1} x_0(r)cas(\frac{2\pi}{N/2}rk) + cos(\frac{2\pi}{N}k)\sum_{r=0}^{\frac{N}{2}-1}x_1(r)cas(\frac{N}{2}rk) \\ + sin(\frac{2\pi}{N}k)\sum_{r=0}^{\frac{N}{2}-1}x_1(r)cas(-\frac{2\pi}{N/2}rk)$$

令$X_{0H}(k) = DHT[x_0(n)]$，$X_{1H}(k) = DHT[x_1(n)]$，可以写成：

$$X_H(k) = X_{0H}(k) + cos(\frac{2\pi}{N}k)X_{1H}(k)+sin(\frac{2\pi}{N}k)X_{1H}(\frac{N}{2}-1)$$

相比于$DIT-DFT$该式中，多了一项$X_{1H}(\frac{N}{2}-k)$。所以可用$X_{0H}(k)$,$X_{1H}(k)$,$X_{0H}(\frac{N}{2}-k)$,$X_{0H}(\frac{N}{2}+k)$四个点同址计算出$X_{H}(k)$,$X_{H}(\frac{N}{2}+k)$,$X_{H}(\frac{N}{2}-k)$,$X_{H}(N-k)$，与$FFT$中的蝶形运算类似，我们叫这种运算“哈特莱蝶形”

令$C(k)=cos(\frac{2\pi}{N}k)$,$S(k)=sin(\frac{2\pi}{N}k)$,当$N \geq 8$时，有

$$X_H(k) = X_{0H}(k) + [C(k)X_{1H}(k)+S(k)X_{1H}(\frac{N}{2}-k)]\\ X_H(\frac{N}{2}+k) = X_{0H}(k) - [C(k)X_{1H}(k)+S(k)X_{1H}(\frac{N}{2}-k)]\\ X_H(\frac{N}{2}-k) = X_{0H}(\frac{N}{2}-k) - [S(k)X_{1H}(k)-C(k)X_{1H}(\frac{N}{2}-k)]\\ X_H(N-k) = X_{0H}(\frac{N}{2}-k) - [S(k)X_{1H}(k)-C(k)X_{1H}(\frac{N}{2}-k)]$$

信号流图

基2 DIT-FHT的运算量(只考虑蝶形变换)

$N = 2^M$
乘法次数：$NM - 3N + 4$
加法次数：$\frac{3}{2}NM - \frac{3}{2}N + 2$

应用

利用$FHT$以及循环卷积定理，与DFT类似，我们就可以通过循换卷积来求出两个实序列的线性卷积。

代码 [UR#34]多项式乘法

#include <cstdio>
#include <cmath>
#include <algorithm>
#define DXT(X,Y) ( XT = (X) , X = ((XT) + (Y)) , Y = ((XT) - (Y)) )
using namespace std;

int N1 , N2 , N4 , L1 , L2 , L3 , L4, n , m , rev[2000001];
double XT , A , E , CC1 , SS1 , T1 , T2 , a[2000001] , b[2000001] , c[2000001];

void FHT( double X[] , int N , int M ) {
    for(int i = 0 ; i < N ; ++i) if(i < rev[i]) swap(X[i] , X[rev[i]]);
    for(int i = 0 ; i < N ; i += 2) DXT(X[i] , X[i+1]);
    N2 = 1; --M;
	if(M <= 0) return;
    while( M-- ) {
        N4 = N2 , N2 = N4 + N4 , N1 = N2 + N2;
        E = 6.283185307179586 / N1;
        for(int j = 0 ; j < N ; j += N1) {
            L2 = j + N2 , L3 = j + N4 , L4 = L2 + N4;
            DXT(X[j] , X[L2]) , DXT(X[L3] , X[L4]) , A = E;
            for(int i = 1 ; i < N4 ; ++i , A += E) {
                L1 = j + i , L2 = j - i + N2;
                L3 = L1 + N2 , L4 = L2 + N2;
				CC1 = cos(A) , SS1 = sin(A);
                T1 = X[L3] * CC1 + X[L4] * SS1;
                T2 = X[L3] * SS1 - X[L4] * CC1;
                XT = X[L1] , X[L1] = XT + T1 , X[L3] = XT - T1;
                XT = X[L2] , X[L2] = XT + T2 , X[L4] = XT - T2;
            }
        }
    }
}

int main() {
    scanf("%d%d" , &n , &m);
    for(int x , i = 0 ; i <= n ; ++i) scanf("%d" , &x) , a[i] = (double)x;
    for(int x , i = 0 ; i <= m ; ++i) scanf("%d" , &x) , b[i] = (double)x;
    int nn = n + m + 1 , L = 0 , tmp;
    for(tmp = 1 ; tmp < nn ; tmp <<= 1, ++L) ; nn = tmp;
    for(int i = 0 ; i < nn ; ++i)
        rev[i] = (rev[i >> 1] >> 1 ) | ((i & 1) << (L - 1));

    FHT(a , nn , L);   FHT(b , nn , L);

	c[0] = b[0] * a[0];
    for(int i = 1; i < nn; ++i) {
		T1 = a[i] , T2 = a[nn - i];
        c[i] = (b[i] * (T1 + T2) + b[nn - i] * (T1 - T2)) * 0.5;
    }
    FHT(c  , nn , L);
	for(int i = 0 ; i <= n + m ; ++i) c[i] /= nn;
    for(int i = 0 ; i <= n + m ; ++i) printf("%d " , (int)(c[i] + 0.5));
    return 0;
}

---

这份代码比我想像的要慢好多，主要原因是循环卷积定理里面运算的增多，以及蝶形变化中较为麻烦的计算四个位置，还有每次使用$cosA$和$sinA$都是很大的开销，曾经尝试使用乘法代替每次计算$cos$和$sin$但是点数$N$很大时，精度问题会被放大的很严重，另一个原因是我不会卡常...但是内存上看起来还是比较优秀的，拜托善于优化的同学优化一下，或者有更好的实现的话，也请告诉我。。。不过看起来还是没人会用啊。。。法法塔的历史地位还是很难被替代的吧。。。改日会补上信号流图。～掰掰