快速傅里叶变换FFT原理及源程序Word下载.docx
《快速傅里叶变换FFT原理及源程序Word下载.docx》由会员分享,可在线阅读,更多相关《快速傅里叶变换FFT原理及源程序Word下载.docx(9页珍藏版)》请在冰点文库上搜索。
![快速傅里叶变换FFT原理及源程序Word下载.docx](https://file1.bingdoc.com/fileroot1/2023-5/5/634dfece-d5c5-4fc7-8f25-ac87aff64b3b/634dfece-d5c5-4fc7-8f25-ac87aff64b3b1.gif)
在输入序列中是按码位倒序排列的,输出序列是按顺序排列;
每级包含个蝶形单元,第级有个群,每个群有个蝶形单元;
每个蝶形单元都包含乘和系数的运算,每个蝶形单元数据的间隔为,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<
n/2;
i++)
{m=i;
j=2*j;
if(j==n)break;
}//得到流程图的共几级
n1=n-1;
for(j=0,i=0;
n1;
{if(i<
j)//如果i<
j,即进行变址
{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;
n;
i+=2)
{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;
i+=n1)//控制同一蝶形结运算,即计算系数相同蝶形结
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<
math.h>
stdio.h>
stdlib.h>
#definepi3.14159265359
{inti,j,k,l,i1,i2,i3,i4,n4,m,n1,n2;
doublea,e,cc,ss,tr,t1,t2;
j)
(j+1))
l++)
i+=n1)
j++)
main()
{FILE*p;
inti,j,n;
doubledt=0.001;
doublex[512];
p=fopen("
d:
\123.c"
"
w"
);
n=512;
{x[i]=sin(200*pi*i*dt);
{fprintf(p,"
%10.7f"
x[i]);
fprintf(p,"
\n"
printf("
fft(x,n);
\nDISCRETEFOURIERTRANSFORM\n"
x[0]);
%10.7f+J%10.7f\n"
x[1],x[n-1]);
for(i=2;
{
%10.7f+J%10.7f"
x[i],x[n-i]);
x[i+1],x[n-i-1]);
x[n/2]);
x[n/2-1],-x[n/2+1]);
x[n/2-i],-x[n/2+i]);
x[n/2-i-1],-x[n/2+i+1]);
将程序运行后所得数据绘制成曲线图(其中FFT变换的数据要先取绝对值后再画图)如下
由上图可知,变换后的图开在频率100Hz处出现一个峰值,这与理论上的结果一致。