快速傅里叶变换FFT原理及源程序之欧阳育创编.docx
《快速傅里叶变换FFT原理及源程序之欧阳育创编.docx》由会员分享，可在线阅读，更多相关《快速傅里叶变换FFT原理及源程序之欧阳育创编.docx（8页珍藏版）》请在冰豆网上搜索。
快速傅里叶变换FFT原理及源程序之欧阳育创编
《测试信号分析及处理》课程作业
时间：
2021.02.04
创作：
欧阳育
快速傅里叶变换
一、程序设计思路
快速傅里叶变换的目的是减少运算量，其用到的方法是分级进行运算。
全部计算分解为
级，其中
；在输入序列
中是按码位倒序排列的，输出序列
是按顺序排列；每级包含
个蝶形单元，第
级有
个群，每个群有
个蝶形单元；每个蝶形单元都包含乘
和
系数的运算，每个蝶形单元数据的间隔为
，i为第i级；同一级中各个群的系数
分布规律完全相同。
将输入序列
按码位倒序排列时，用到的是倒序算法——雷德算法。
自然序排列的二进制数，其下面一个数总比上面的数大1，而倒序二进制数的下面一个数是上面一个数在最高位加1并由高位向低位仅为而得到的。
若已知某数的倒序数是
，求下一个倒序数，应先判断
的最高位是否为0，与
进行比较即可得到结果。
如果
，说明最高位为0，应把其变成1，即
，这样就得到倒序数了。
如果
，即
的最高位为1，将最高位化为0，即
，再判断次高位；与
进行比较，若为0，将其变位1，即
，即得到倒序数，如果次高位为1，将其化为0，再判断下一位……即从高位到低位依次判断其是否为1，为1将其变位0，若这一位为0，将其变位1，即可得到倒序数。
若倒序数小于顺序数，进行换位，否则不变，防治重复交换，变回原数。
注：
因为0的倒序数为0，所以可从1开始进行求解。
二、程序设计框图
（1）倒序算法——雷德算法流程图
（2）FFT算法流程
三、FFT源程序
voidfft（x,n）
intn;
doublex[];
{inti,j,k,l,m,n1,n2;
doublec,c1,e,s,s1,t,tr;
for（j=1,i=1;i{m=i;
j=2*j;
if（j==n）break;
}//得到流程图的共几级
n1=n-1;
for（j=0,i=0;i{if（i{tr=x[j];
x[j]=x[i];
x[i]=tr;
}
k=n/2;//求j的下一个倒位序
while（k<（j+1））//如果k<（j+1）,表示j的最高位为1
{j=j-k;//把最高位变成0
k=k/2;//k/2，比较次高位，依次类推，逐个比较，直到某个位为0
}
j=j+k;//把0改为1
}
for（i=0;i{tr=x[i];
x[i]=tr+x[i+1];
x[i+1]=tr-x[i+1];
}
n2=1;
for（l=1;l<=m;l++）//控制蝶形结级数
{n4=n2;
n2=2*n4;
n1=2*n2;
e=6.28318530718/n1;
for（i=0;i{tr=x[i];
x[i]=tr+x[i+n2];
x[i+n2]=tr-x[i+n2];
x[i+n2+n4]=-x[i+n2+n4];
a=e;
for（j=2;j<=（n4-1）;j++）//控制计算不同种蝶形结，即计算系数不同的蝶形结
{i1=i+j;
i2=i-j+n2;
i3=i+j+n2;
i4=i-j+n1;
cc=cos（a）;
ss=sin（a）;
a=a+e;
t1=cc*x[i3]+ss*x[i4];
t2=ss*x[i3]-cc*x[i4];
x[i4]=x[i2]-t2;
x[i3]=-x[i2]-t2;
x[i2]=x[i1]-t1;
x[i1]=x[i1]+t1;
}
}
}
}
四、计算实例及运行结果
设输入序列
为
其离散傅里叶变换为
这里
。
选n=512,计算离散傅里叶变换
。
所用软件为Turboc2.0，操作界面如图1所示
图1Turboc2.0操作界面
程序运行结束后的界面如图2所示
图2程序运行后的界面
例子的具体程序如下：
#include
#include
#include
#definepi3.14159265359
voidfft（x,n）
intn;
doublex[];
{inti,j,k,l,i1,i2,i3,i4,n4,m,n1,n2;
doublea,e,cc,ss,tr,t1,t2;
for（j=1,i=1;i{m=i;
j=2*j;
if（j==n）break;
}
n1=n-1;
for（j=0,i=0;i{if（i{tr=x[j];
x[j]=x[i];
x[i]=tr;
}
k=n/2;
while（k<（j+1））
{j=j-k;
k=k/2;
}
j=j+k;
}
for（i=0;i{tr=x[i];
x[i]=tr+x[i+1];
x[i+1]=tr-x[i+1];
}
n2=1;
for（l=1;l<=m;l++）
{n4=n2;
n2=2*n4;
n1=2*n2;
e=6.28318530718/n1;
for（i=0;i{tr=x[i];
x[i]=tr+x[i+n2];
x[i+n2]=tr-x[i+n2];
x[i+n2+n4]=-x[i+n2+n4];
a=e;
for（j=2;j<=（n4-1）;j++）
{i1=i+j;
i2=i-j+n2;
i3=i+j+n2;
i4=i-j+n1;
cc=cos（a）;
ss=sin（a）;
a=a+e;
t1=cc*x[i3]+ss*x[i4];
t2=ss*x[i3]-cc*x[i4];
x[i4]=x[i2]-t2;
x[i3]=-x[i2]-t2;
x[i2]=x[i1]-t1;
x[i1]=x[i1]+t1;
}
}
}
}
main（）
{FILE*p;
inti,j,n;
doubledt=0.001;
doublex[512];
p=fopen（"d:
\123.c","w"）;
n=512;
for（i=0;i{x[i]=sin（200*pi*i*dt）;
}
for（i=0;i{fprintf（p,"%10.7f",x[i]）;
fprintf（p,"\n"）;
printf（"%10.7f",x[i]）;
printf（"\n"）;
}
fft（x,n）;
fprintf（p,"\nDISCRETEFOURIERTRANSFORM\n"）;
printf（"\nDISCRETEFOURIERTRANSFORM\n"）;
fprintf（p,"%10.7f",x[0]）;
printf（"%10.7f",x[0]）;
fprintf（p,"%10.7f+J%10.7f\n",x[1],x[n-1]）;
printf（"%10.7f+J%10.7f\n",x[1],x[n-1]）;
for（i=2;i{
fprintf（p,"%10.7f+J%10.7f",x[i],x[n-i]）;
fprintf（p,"%10.7f+J%10.7f",x[i+1],x[n-i-1]）;
fprintf（p,"\n"）;
printf（"%10.7f+J%10.7f",x[i],x[n-i]）;
printf（"%10.7f+J%10.7f",x[i+1],x[n-i-1]）;
printf（"\n"）;
}
fprintf（p,"%10.7f",x[n/2]）;
printf（"%10.7f",x[n/2]）;
fprintf（p,"%10.7f+J%10.7f\n",x[n/2-1],-x[n/2+1]）;
for（i=2;i{
fprintf（p,"%10.7f+J%10.7f",x[n/2-i],-x[n/2+i]）;
fprintf（p,"%10.7f+J%10.7f",x[n/2-i-1],-x[n/2+i+1]）;
fprintf（p,"\n"）;
printf（"%10.7f+J%10.7f",x[n/2-i],-x[n/2+i]）;
printf（"%10.7f+J%10.7f",x[n/2-i-1],-x[n/2+i+1]）;
printf（"\n"）;
}
}
将程序运行后所得数据绘制成曲线图（其中FFT变换的数据要先取绝对值后再画图）如下
由上图可知，变换后的图开在频率100Hz处出现一个峰值，这与理论上的结果一致。
时间：
2021.02.04
创作：
欧阳育
升级会员 