1、FFT
本文主要依据前面所写的文章请戳进行规律总结,便于后续的程序实现。这里仍然以8点fft进行分析,具体示意图如下:
由上图的输入分析:
输入为 x ( 0 ) , x ( 4 ) , x ( 2 ) , x ( 6 ) , x ( 1 ) , x ( 5 ) , x ( 3 ) , x ( 7 ) x(0),x(4),x(2),x(6),x(1),x(5),x(3),x(7) x(0),x(4),x(2),x(6),x(1),x(5),x(3),x(7),这些输入数据的序号并不是按照顺序的;
将 x ( 0 ) − x ( 7 ) x(0)-x(7) x(0)−x(7)用二进制表示:
x ( 000 ) , x ( 001 ) , x ( 010 ) , x ( 011 ) , x ( 100 ) , x ( 101 ) , x ( 110 ) , x ( 111 ) x(000),x(001),x(010),x(011),x(100),x(101),x(110),x(111) x(000),x(001),x(010),x(011),x(100),x(101),x(110),x(111),
将二进制进行翻转,得到:
x ( 000 ) , x ( 100 ) , x ( 010 ) , x ( 110 ) , x ( 001 ) , x ( 101 ) , x ( 011 ) , x ( 111 ) x(000),x(100),x(010),x(110),x(001),x(101),x(011),x(111) x(000),x(100),x(010),x(110),x(001),x(101),x(011),x(111)
对应的十进制序号为:
x ( 0 ) x ( 4 ) ∣ x ( 2 ) x ( 6 ) ∣ x ( 1 ) x ( 5 ) ∣ x ( 3 ) x ( 7 ) x(0) x(4) | x(2) x(6) | x(1) x(5) | x(3) x(7) x(0)x(4)∣x(2)x(6)∣x(1)x(5)∣x(3)x(7)
这个正好是上面按照奇偶分开得到的输入序列,具体程序实现如下。
void bitReverse(float* xreal, float* ximag, int n)
{
int i = 0, j = 0, a = 0, b = 0, p = 0;
for (i = 1; i < n; i *= 2)
p++;
for (i = 0; i < n; i++)
{
a = i;
b = 0;
for (j = 0; j < p; j++)
{
b = (b << 1) + (a & 1);
a >>= 1;
}
if (b > i)
{
swap(&xreal[i], &xreal[b]);
swap(&ximag[i], &ximag[b]);
}
printf("b = %d\n", b);
}
}
说明:
(1)m代表某一级(N = 8时,一共可以分为m = 0,1,2三级),M表示一共有多少级
(2)旋转因子记为: W N p W_{N}^{p} WNp,p是旋转因子指数
(3)旋转因子的增量记为:k,比如在m级最后一级增量k总为1
(4)B表示每个蝶形运算的输入数据的间隔
由m-1级到m级蝶形运算可表示为:
A ( j ) m = A ( j ) m − 1 + A ( j + B ) m − 1 ∗ W N p A(j)_{m} = A(j)_{m - 1}+A(j+B)_{m - 1} * W_{N}^{p} A(j)m=A(j)m−1+A(j+B)m−1∗WNp
A ( j + B ) m = A ( j ) m − 1 − A ( j + B ) m − 1 ∗ W N p A(j+ B)_{m} = A(j)_{m - 1} - A(j + B)_{m - 1} * W_{N}^{p} A(j+B)m=A(j)m−1−A(j+B)m−1∗WNp
对于c实现,由于c没法直接进行复数运算,可以进一步推导如下:
由图一分析可知
(1)每个蝶形运算中两个输入数据的间隔B = 2 m 2^{m} 2m
(2)有 2 m 2^{m} 2m个旋转因子
(3)旋转因子 W W W增量k = 2 M − m 2^{M - m} 2M−m
(4)同一W的间隔为: 2 m + 1 2^{m + 1} 2m+1
(5)蝶形运算的种类数目 2 m 2^{m} 2m等于间隔B
(6)同种蝶形运算次数k = 2 M − m 2^{M - m} 2M−m
下面就可以进行代码实现了,具体实现主要分为三个循环:(1)fft的级数 (2)每一级dft的次数 (3)每一个蝶形运算;
void fft(float* xreal, float* ximag, int n, int M)
{
bitReverse(xreal, ximag, n);
int m = 0, B = 1, j = 0, k = 0, P = 0, i = 0, r = 0;
float Tr = 0, Ti = 0, theta;
theta = 2.0 * PI / n;
for (m = 0; m <= M; m++) //一共有多少级
{
B = 1;
B = (int)pow(2, m);
printf("---------------m=%d---------------\n",m);
for (j = 0; j <= B - 1; j++) //一共有多少个dft
{
k = (int)pow(2, M - m);
P = 1;
P = j * k;
for (i = 0; i <= k - 1; i++) //每个dft的运算
{
r = j + 2 * B * i;
printf("(%d,%d)\n", r, r + B);
Tr = xreal[r + B] * cos(P * theta) + ximag[r + B] * sin(P * theta);
Ti = xreal[r + B] * cos(P * theta) - ximag[r + B] * sin(P * theta);
xreal[r + B] = xreal[r] - Tr;
ximag[r + B] = ximag[r] - Ti;
xreal[r] = xreal[r] + Tr;
ximag[r] = ximag[r] + Ti;
}
}
}
}
2、IFFT
知道了fft的运算方式,求IFFT,我们可以从一下两个方面来进行求解:
ifft的实现:
void ifft(float *xreal, float *ximag, int n, int M)
{
for (int k = 0; k < n; k++)
{
ximag[k] = -ximag[k];
}
fft(xreal, ximag, n, M);
for (int k = 0; k < n; k++)
{
xreal[k] = xreal[k] / n;
ximag[k] = -ximag[k] / n;
}
}
水平有限,若有不当之处请指教。