FFT算法设计含程序设计Word文档下载推荐.docx
- 文档编号:16259099
- 上传时间:2022-11-22
- 格式:DOCX
- 页数:16
- 大小:541.44KB
FFT算法设计含程序设计Word文档下载推荐.docx
《FFT算法设计含程序设计Word文档下载推荐.docx》由会员分享,可在线阅读,更多相关《FFT算法设计含程序设计Word文档下载推荐.docx(16页珍藏版)》请在冰豆网上搜索。
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⑶
A⑷
A⑸
A(7)
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;
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
000
001
100
2
3
010
011
110
6
5
7
100^□1011
111
(001
101
倒序规律
对于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<
N/4),则直接加1(J<
=J+N/4),否则先将次高位变成0(J<
=J-N/4),再判断下一位,依次类推,直到完成高位加1,适2(1+1二10b,向右进位的运算。
LH<
=N/2
J<
=LH
I可以从1变换到(N/2-1),
Nk=N-2
三I
<
n>
J[即后半部不用考虑
」丄如果该高位为1,则先减去N/2或N/4、N/8...再判断下个次高位
//输入序列倒序软件程序
j=N/2;
//第0个数(二进制数都为0)和最后一个第N4个数(二进制数都为丄)不需倒序for(i=1;
i<
N-2;
i++)
temp=dataR[i];
dataR[i]=dataR[j];
dataR[j]=temp;
}
k=N/2;
while(l)
输入序列倒序的算法设计方法二
十进制数X
O1O
W厶2
W厶V
ni1
11n
KJXX
XXKJ
w
001log
Oil
从表格可以看出,所谓倒序只是把数组下标的……
最高位与最低位互换次高位与次低位互换
方法二软件分析:
已一个字节(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;
//取次低位
0x40;
//取次高位
if(F0)j=j|0x40;
//最次位换到次高位if(Fl)j=j|0x02;
//最次位换到次低位
0x04;
Fl=i&
0x20;
if(FO)j=j|0x20;
if(Fl)j=j|0x04;
0x08;
0xl0;
if(FO)j=j|OxlO;
if(Fl)j=j|0x08;
//前半部与后半部交换,相等时无需交换
A[i]OAU]
■
9
只需单层循环即可,共需要循环(128-2)次
算法改进一:
前面的算法可以进一步优化为:
for(i=0;
i++)
for(k=0;
k<
4;
k++)
(0x01<
k);
(0x80>
>
if(FO)j=j|(0x80»
//最低位换到最高位iff(Fl)j=j|(OxOl«
//最高位换到最低位
if(i<
j)//前半部与后半部交换,相等时无需交换
A[i]OA[j];
这个算法只是针对8位的,如果是任意N位的应该如何做?
这里的N必须满足N=2M
算法改进二:
针对任意N=2M的情况
(N/2);
i++)//或#for(i=0;
((N-l)/2);
i++)
j=o;
(M/2);
k++)//当N=256时,M=8
m=0x01<
k;
n=0x80>
FO=i&
m;
Fl=i&
n;
if(FO)j=j|n;
if(Fl)j=j|m;
FFT软件示例
#include<
math.h>
#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;
N;
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
//输入序列倒序
//第0个数(二进制数都为0)和最后一个第N7个数(二进制数都为丄)不需倒序for(i=1;
i++){
if(i<
J)<
dataRU]=temp;
//因为波形虚数部分都为6所以不用交换
//temp=datal[i];
//datal[i]=datal[j];
//datalO]=temp;
while(l){
(j<
k)<
j=j+k;
break;
}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;
=B-l;
J++)//J=-1
//计算P=J*2^(M-L)
p=i=M-L;
while(i)<
p=p*2;
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;
N/2;
{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 算法 设计 程序设计
![提示](https://static.bdocx.com/images/bang_tan.gif)