标签：name 正交 isp 没有 西安电子 数字信号 加速 const std
最近复习\(DSP\)的时候，发现了一个号称专门针对离散实序列的变换，经分析总运算量为普通\(FFT\)的几乎一半，而且完全没有复数。这么强的吗？于是花了一个下午，去学习了一下。中文资料少的可怜。。。于是去读书馆翻了几乎所有的\(dsp\)课本。。。发现了这本书 西安电子科技大学出版社《数字信号处理》第二版！竟然花了一节在讲\(DHT\)和\(FHT\)！竟然还有附录FORTRAN代码（虽然一堆错）！这里把这个东西稍微普及一下。。。不过估计没人用的上
对于序列\(x(n)\)其\(N\)点\(DFT\)具有简单的共轭对称性，即\[X(N - k) = X^*(k), k = 0, 1 ... ,N-1 \]
所以只要计算\(X(k)\)的前\(N/2\)个值，则后\(N/2\)可以通过\((1)\)求得，\(X(k)\)的\(N/2\)个复数正好对应\(N\)个实数数据(\(N\)点实序列\(DFT\)，也可以通过把\(N\)个实数，压成\(N/2\)个复数，再利用共轭对称性求解，这里不做详解)。由此可见，一个\(N\)点实序列的\(DFT\)，完全可由\(N\)个实数数据确定。由此，我们引出一种直接对于实序列进行实数域变换的离散哈特莱变换(\(DHT\))。\(DHT\)将\(x(n)\)看成实数数据，而不是像\(DFT\)一样看作虚部为\(0\)的复数，因此节省了一半的空间，运算效率也提升了近一倍，且\(DFT\)与\(DHT\)之间存在简单的关系，容易实现相互转换。
设\(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. \]
用\(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\)(隐含周期性)。
若\(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\)的运算更加复杂，这一点把这个算法的在竞赛中的实用性大大降低了。。。毕竟我就是为了加速卷积，不过他空间上的优势还是很名明显的。
与\(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)]
\]
\(N = 2^M\)
乘法次数：\(NM - 3N + 4\)
加法次数：\(A_H = \frac{3}{2}NM - \frac{3}{2}N + 2\)
利用\(FHT\)以及循环卷积定理，与DFT类似，我们就可以通过循换卷积来求出两个实序列的线性卷积。
#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的历史地位还是很难被替代的吧。。。改日会补上信号流图。掰掰
标签：name 正交 isp 没有 西安电子 数字信号 加速 const std
原文地址：https://www.cnblogs.com/RRRR-wys/p/10090007.html
