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;
}
}
水平有限,若有不当之处请指教。