離散哈特萊變換(DHT)及快速哈特萊變換(FHT)學習

說在前邊

最近複習\(DSP\)的時候，發現了一個號稱專門針對離散實序列的變換，經分析總運算量為普通\(FFT\)的幾乎一半，而且完全沒有複數。這麼強的嗎？於是花了一個下午，去學習了一下。中文資料少的可憐。。。於是去讀書館翻了幾乎所有的\(dsp\)課本。。。發現了這本書 西安電子科技大學出版社《數字訊號處理》第二版！竟然花了一節在講\(DHT\)和\(FHT\)！竟然還有附錄FORTRAN程式碼（雖然一堆錯）！這裡把這個東西稍微普及一下。。。不過估計沒人用的上

離散哈特萊變換(DHT)

引入

對於序列\(x(n)\)其\(N\)

預備知識

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)\)

\[ \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\)的轉換容易實現

DFT的性質

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)\),\(H_{0H}(\frac{N}{2}-k)\),\(H_{0H}(\frac{N}{2}+k)\)四個點同址計算出\(X_{H}(k)\),\(X_{H}(\frac{N}{2}+k)\),\(H_{H}(\frac{N}{2}-k)\),\(H_{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\)
加法次數：\(A_H = \frac{3}{2}NM - \frac{3}{2}N + 2\)

應用

利用\(FHT\)以及迴圈卷積定理，與DFT類似，我們就可以通過循換卷積來求出兩個實序列的線性卷積。

程式碼 [UR#34]多項式乘法

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

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

void FHT( double X[] , int N , int M , int f ) {
    for(int i = 0 ; i < N ; ++i) if(i < rev[i]) {
        XT = X[i] , X[i] = X[rev[i]] , X[rev[i]] = XT;
    }
    for(int i = 0 ; i < N ; i += 2) DXT(X[i] , X[i+1]);
    N2 = 1;
    for(int k = 1 ; k < M ; ++k) {
        N4 = N2 , N2 = N4 + N4 , N1 = N2 + N2;
        E = PI * 2.0 / 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;
            }
        }
    }
    if(f == -1) {
        for(int i = 0 ; i < N ; ++i) X[i] /= N;
    }
}

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;
    int tmp; L = 0;
    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 , 1), FHT(b , nn , L , 1);
    for(int i = 0; i < nn; ++i) {
        e = i > 0 ? (a[i] + a[(nn - i)]) * 0.5 : a[0];
        o = i > 0 ? (a[i] - a[(nn - i)]) * 0.5 : 0.0;
        c[i] = b[i] * e + b[i == 0 ? 0 : (nn - i)] * o;
    }
    FHT(c  , nn , L , -1);
    for(int i = 0; i <= n+m; ++i) printf("%d ",(int)(c[i]+0.5));
    return 0;
}

這份程式碼比我想像的要慢好多，主要原因還是迴圈卷積定理裡面運算的增多，另一個原因是我不會卡常...但是記憶體上看起來還是比較優秀的，拜託善於優化的同學優化一下，或者有更好的實現的話，也請告訴我。。。不過看起來還是沒人會用啊。。。fft的歷史地位還是很難被替代的吧。。。改日會補上訊號流圖。掰掰