FFT算法设计含程序设计.docx
- 文档编号:3353601
- 上传时间:2022-11-22
- 格式:DOCX
- 页数:16
- 大小:541.44KB
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(J LH<=N/2 J<=LH I可以从1变换到(N/2-1), Nk=N-2 三I 」丄如果该高位为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(i A[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(i A[i]OA[j]; FFT软件示例 #include #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;i if(i temp=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=J*2^(M-L) 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
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- FFT 算法 设计 程序设计