FFT算法设计含程序设计.docx
《FFT算法设计含程序设计.docx》由会员分享，可在线阅读，更多相关《FFT算法设计含程序设计.docx（16页珍藏版）》请在冰豆网上搜索。
FFT算法设计含程序设计
第6讲FFT算法设计
傅立叶变换将信号从时域转换为频域，可以进行模拟信号的频率分析
离散傅立叶变换（DFT）将信号从频域转换为数字（频）域，可以进行数字信号（模拟信号数字化）的频率分析
为了实现DFT在计算机上的快速实现，提出了快
速离散傅立叶变换（FFT）
如何有傅氏变换.>DFT->FFT?
ej（,）=cosco±jsin
欧拉公式：
尸0）=匚f^e-jwtdt=[f（t）e~^dt=f（t）e~^dt
-0kn
N7.j二
=>X（k）=A（/z）e“仏=0丄2,N
fi=0
令必=「不F称为旋转因子
N_\
X伙）=工x（n）W^n,k=0,1,2,……“—1
上式中，k对应数字域，n对应时域
另有推导时需用到的公式：
1）W'^+IN=e小'
1N为I个周期
—j——m
=eN=W；1
7、-7—*（-w）
Z）乂丁=en
=W^m
N2兀
3）%2一
I对称性
周期性
j聲e，N・m为加上一个周期
|周期性僅韦
—,N、
=-®
4）叽一F
-j—m一-J
eN*eN2=€
，其中e~j7T—cos（tt）—Jsin（^）=-1
・2兀
-J——*.N_
一"n
可约哲
推导分析
若序列x（n）的长度为N,且满足N=2M,（M为自然数）按n的奇偶性把x（n）分解为两个N/2的子序列：
x1（r）=x（2r）zr=0zlz..vN/2-1x2（r）=x（2r+l）zr=Oflz..vN/2-1
则x（n）的DFT为：
X（R）=工.心）歐”+工
，匸偶数，匸奇数g
N/2-1N/2-I
=^x（2r）W~kr+为班2厂+l）W?
m）
r=（）r=（）
N/2-1N/2-1
=£石（厂）£-勺什）吒"
r=0r=0
r=0
其中XJk）和X2（k）均以N/2为周期
NN梯N
X（k+于）=X|伙+亍）+叫2X?
伙+刁）
同理，可推出：
x"）=X?
伙）+W爲总4仏）
n，k=Ozlz...,N/4-1
X^+-）=X.（k）-W^,2X4（k）
4
X?
伙）—X5伙）十w爲2乂6伙）
x十牛=X5（k）-W爲2X9k=0几…,N/4-1
分到最后，k=0,进行蝶形运算的两个输入就是最初输入序列x（n）的其中两个。
X（0）x（l）X⑵X⑶X（4）X⑸X（6）X（7）
N=8点FFT运算图示
x（0）
"4）3
M2）
%（6）f竺）
x（l）M也
x（5）朋
x（3）3也
x（7）皿辺
A（0）A（0）
A（0）
A（l）
A
（2）
A
（2）
A⑶
A⑷
A⑸
A⑸
A（7）
A（l）
A（6）
妙x（o）AOlx（l）
X⑵迥X⑶A（4x⑷^X（5）^X（6）&Zlx（7）
N=16点FFT运算图示
蝶形运算规律
设序列x（n）已经经过时域抽选（倒序）后，存入数组X中。
如果蝶形运算的两个输入相距B个点，用原位计算（即计算结果还放在数组的原来位置），则蝶形运算可表示成如下形式：
X伙）=X（伙）十呎X2伙）
X（Z普）=X’（約一比兀（幻
X,（丿）＜=X-（丿）+X」丿+B）Wf；⑴
X/（y+〃）v=X/」gz）—X,+（并
其中：
p=J*2ML;J=0,lf...f2L1-lL=1Z2,...,M
下标L表示第L级运算，XL（J）则表示第L级运算后数组元素XQ）的值。
如果要用实数运算完成上述蝶形运算，可由下面的算法进行:
设：
T=XZ=TR+jTl（3）下标R表示实部
X/1（J）=X/?
（/）+jX,J）（4）下标1表示虚部X）Q）代表上一次的实数值
X-（八3）=（八8）+JX/（J+B）（5）X_G/+3）W/=[XQ+3）fX/G/+3）]W/
=lXK（J+B）+jX/（J+B）|（cos等p-./sin等/?
）
=[XR（J+B）cos—p+X,（J+B）sin—p]+j[X,（J+B）cos"pXN（J+B）sin"p\NNNN
Tr=Xr（J+B）cos#p+x；（丿+B）sin千p7；=X/（丿+B）cos#p一Xr（丿+5）sin—p
XL（J）=XR（J）+jX/（J）
=XL_SJ）+TR+jT{由⑴（6）式推导得出
=X'RU）+jX\m+TR+丿刀由⑷式推导得出
X.（丿+3）=X-（J）一（7；+JTf）由
（2）（6）式推导得出
=X、RQ）+jX,J）-（jR+j"）由（4）式得出
=>厂应心乂斥⑴+几
"[X/（J）=X；（J）+T/
Xr（J+B）=Xk（J）-Tr
X/（J+B）=X/（J）-T/
输入序列倒序
公式（7）（8）（9）主要用于FFT的软件编程
Tr=X,八〃）cos#“+X/（丿+B）sin#
T,=X/（丿+B）cos-^p-XN（J+B）sin二p
Xr（J）=Xr（J）+Tr
X/（J）=X；（J）+T/
Xr（J+B）=Xr（J）—Tr
B<=2A（L-1）
P=J*2A（M-L）
X伙）v=X（A）+X（Ar+B）W^X（k+B）<=Xg-X（k+B）W；
Xf（J+B）=Xl（J）-T/
公式中的J就是流程图中公式的变量k,流程图中：
N表示阶数，
M表示总级数，
L表示当前级数，
B表示每个蝶形的两个输入数据的间隔,
P表示旋转因子指数可以看出，流程图总共有3个循环外循环：
次数为级数L的变换范围中循环：
为根据当前L求出各个不同的P，循环次数为P的个数2-丄内循环：
为每级中每个P对应的蝶形运算个数,循环次数为2M』
内循环中k值每次变换范围（增量）为
这是同一级中每个相同的P对应不同蝶形运算的间隔。
看图推导软件编程规则：
方法
1.当N=8时，第L级共有2-丄个不同的旋转因子。
（P称为旋转因子指数）k=2（k为p的增量）
因为N=2M,所以有L=1Z2Z...ZM,即L的最大值为M当L"时,p=0；当1_=2时/p=°,2;
当L=3时‘pnOhaQ;k=l
当N=16ffif
当L=J■时.p=O;
当L=2时,p=O,4;k=4
当L=3时,p=0,2,4,6;k=2
当L=4时/p=0flz2F3f4f5F6z7;k=l
2・14/卩=/尸*（j=6丄‘2,3/…）
（归纳得出：
N/k=2L）
=（L=lz2z3f…表当前级数）
=wj
（M表总级数）
3■每个p对应每级中的运算个数为:
L=4=M
L=3
L=1L=2
8个块
pi=8
4个蝶形块有2个蝶形块
pi=4
pi=2
有1个蝶形块，pi=C
pi定义为p的增竄
M-4=2^-Af=2°
当L=M时，pi=l=2°当厶=M-1H寸，pi=2=2*当L=M2时，pi=4=22
当厶=Mr3=l时，pi=8=2
'当乙=4H寸，pi=2当L=3H寸，pi=2M3当L=2H寸，pi=2M2当厶=l时，pi=2""
=>pi=2M_L
令p=J*pi=J*2ML（其中J=0,1,2,3,…），这样写是为了利于软件
的循环编程。
此时只要求出每级J的个数J”⑹即可
%=4时，JTolaI=8=2
J的总个数J“讪为2",每一级P的总个数也为2"
外循环次数为级数L
中循环为根据当前L求出各个不同的p,循环次数为p的个数2-丄
内循环为每级中每个p对应的蝶形运算个数（记为⑹），循环次数为2"丄
当乙=1时，VTixal=8=23
当乙=2时，Vt如=4=2?
当乙=3时，^.^=2=2*
当无=4时，J侖=1=2°
每个蝶形的两个输入数据间隔（记为INd）:
当厶=1吋，LV^=1=2°'当£=2H寸，IN（I=2=2'当厶=3时，Wd=4=22当厶=4时，Wtl=8=23
同一级中每个相同的P对应蝶形运算的间隔（记为VQ：
=>Vd=2L
当L=W寸,Vz=2=2,当L=2B寸,Vz=4=2?
当L=3H寸，Vz=8=23当厶=4（甘，V〃=16=24
可以看出，为了利于编程，所有的公式推导都往L上靠
输入序列倒序的算法设计
顺序
倒序
十进制数I
二进制数
二进制数
十进制数J
0
000
000
0
1
001
100
4
2
3
010
011
010
110
2
6
4
5
6
7
100^□1011
110
111
（001
101
011
111
1
5
3
7
倒序规律
对于N=2M,M位二进制数最高位的权值为N/2,且从左向右二进制位的权值依次为N/4,N/8,...,2,1。
因此，最高位加1相当于十进制运算J+N/2o如果最高位是0（JvN/2）,则直接由J+N/2得下一个倒序值；如果最高位是1（JNN/2）,则要将最高位变成O（Jv二J・N/2）,次高位加1（J+N/4）。
但次高位加1时同样要判断0、1值，如果为0（JLH<=N/2
J<=LH
I可以从1变换到（N/2-1）,
Nk=N-2
三I
J［即后半部不用考虑
」丄如果该高位为1,则先减去N/2或N/4、N/8...再判断下个次高位
//输入序列倒序软件程序
j=N/2;
//第0个数（二进制数都为0）和最后一个第N4个数（二进制数都为丄）不需倒序for（i=1;i<
temp=dataR[i];dataR[i]=dataR[j];dataR[j]=temp;
}
k=N/2;
while（l）
输入序列倒序的算法设计方法二
十进制数X
二进制数
二进制数
十进制数J
000
000
0
1
001
100
4
2
O1O
O1O
2
W厶2
W厶V
ni1
11n
KJXX
XXKJ
w
6
7
100
101
110
111
001log
Oil
111
1
5
3
7
从表格可以看出，所谓倒序只是把数组下标的……
最高位与最低位互换次高位与次低位互换
方法二软件分析：
已一个字节（N=256）的倒序为例
A[0]zA[l]zA[255]
（下标从0000_0000变化到1111_1111）
/*定义两个标志位FO,F1*/for（i=0;i<（255/2）;i++）//除2是因为只要判断前半部{j=o；
F0=i&0x01;//取最低位
Fl=i&0x80;//取最高位
if（FO）j=j|0x80;//最低位换到最高位if（Fl）j=j|0x01;〃最高位换到最低位
F0=i&0x02;//取次低位
Fl=i&0x40;//取次高位
if（F0）j=j|0x40;//最次位换到次高位if（Fl）j=j|0x02;//最次位换到次低位
F0=i&0x04;Fl=i&0x20;if（FO）j=j|0x20;if（Fl）j=j|0x04;
F0=i&0x08;
Fl=i&0xl0;if（FO）j=j|OxlO;if（Fl）j=j|0x08;
//前半部与后半部交换，相等时无需交换
A[i]OAU]
■
9
只需单层循环即可，共需要循环（128-2）次
算法改进一：
前面的算法可以进一步优化为:
for（i=0;i<（255/2）;i++）
for（k=0;k<4;k++）
F0=i&（0x01<Fl=i&（0x80>>k）;//取最高位
if（FO）j=j|（0x80»k）;//最低位换到最高位iff（Fl）j=j|（OxOl«k）;//最高位换到最低位
if（iA[i]OA[j]；
这个算法只是针对8位的，如果是任意N位的应该如何做?
这里的N必须满足N=2M
算法改进二：
针对任意N=2M的情况
for（i=0;i<（N/2）;i++）//或#for（i=0;i<（（N-l）/2）;i++）
for（k=0;k<（M/2）;k++）//当N=256时，M=8
m=0x01<n=0x80>>k;FO=i&m;Fl=i&n;if（FO）j=j|n;
if（Fl）j=j|m;
if（iA[i]OA[j]；
FFT软件示例
#include#definePI3.1415926#defineN128#defineM7
#defineAO255.0
voidmain（void）
intizjzk;
intp,J,L,B;
floatSinIn[N];
floatdataR[N],dataI[N];floatdataA[N];
floatTr^Thtemp;
//输入波形
for（i=0;i<
Sinln[i]=AO*（sin（2*PI*i/25）+sin（2*PI*i*0.4））;dataR[i]=Sinln[i];//输入波形的实数部分
datal[i]=0;//输入波形的虚数部分
dataA[i]=0;//输出波形的幅度谱为0
//输入序列倒序
j=N/2;
//第0个数（二进制数都为0）和最后一个第N7个数（二进制数都为丄）不需倒序for（i=1;iif（itemp=dataR[i];dataR[i]=dataR[j];dataRU]=temp;
//因为波形虚数部分都为6所以不用交换
//temp=datal[i];
//datal[i]=datal[j];//datalO]=temp;
k=N/2;
while（l）{
i<（j}else{j=j-k;
k=k/2;
>
//进行FFT
//Xr[J]=Xrf（J）+Tr
//Xi[J]=X「Q）+Ti//Xr[J+B]=Xrf（J）-Tr//Xi[J+B]=X『Q）・Ti
//（其中X严为上一级的Xivxr为上一级的Xi）
//其中T「=Xr,（J+B）cos（2.0*PI*p/N）+Xi,（J+B）sin（2.0*PI*p/N）
//Ti=Xi,（J+B）cos（2.0*PI*p/N）・Xrf（J+B）sin（2.0*PI*p/N）
for（L=l;L<=M;L++）//FFT蝶形级数L从丄一M
{
//计算每个蝶形的两个输入数据相距B=2^（L-1）;
B=1;
i=L-l;while（i）{
B=B*2;
}//或采用运算，即B=2^（L-1）;B=（int）pow（2,L-l）；
//第L级蝶形有pow（2zL-l）r即2的-丄次方个蝶形运算和pow（2zL-l）个旋转因子pfor（J=0;J<=B-l;J++）//J=-1
p=i=M-L;
while（i）p=J*p;
//第L级蝶形中同一个旋转因子包含多少个蝶形运算
//每个蝶形的两个输入数据相距B=2^（L-1）个点//同一旋转因子对应着间隔为2八L点的2八（M-L）个蝶形（2八1_=2*B）for（k=J;k<=N-l;k=k+2*B）{
Tr=dataR[k];
Ti=datal[k];
temp=dataR[k+B];
dataR[k]=dataR[k]+dataR[k+B]*cos（2.0*PI*p/N）
+dataI[k+B]*sin（2.0*PI*p/N）;
datal[k]=datal[k]+dataI[k+B]*cos（2.0*PI*p/N）
・dataR[k+B]*sin（2.0*PI*p/N）;
dataR[k+B]=Tr・dataR[k+B]*cos（2.0*PI*p/N）
・dataI[k+B]"hi（2・0*PI*p/N）;
dataI[k+B]=Ti・dataI[k+B]*cos（2.0*PI*p/N）
+temp*sin（2.0*PI*p/N）;
//求出幅度频率谱
//因为从0-N是0-2PI范围,所以只要求出0-N/2即可
//幅频对应的位置可由丄28*（1/2PI）=（128/25）*（1/x）求出
]]■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■
//IIII
//（采样频率）丄屏总个数周期（输入频率）丄屏总个数输入频率周期
〃得出输入频率x=2PI/25
〃对应的幅频值的波形位置=128/25=5.12（因为2PI对应点的位置为丄28,PI对应点的位置为64）
for（i=0;i{dataA[i]=sqrt（dataR[i]*dataR[i]+dataI[i]*dataI[i]）;
}
〃如果要算相位，则e[i]=arctan（datal[i]/dataR[i]）
//频谱的平方称为功率谱，记为：
P[i]=dataR[i]*dataR[i]+dataI[i]*dataI[i]while
（1）;//breakpoint
升级会员 