又苦看了两天书,倒是出来了<br />作了三个序列的实验,均作的8点的:<br />x1=2*i+j0 // 0<=i<=7<br />x2=3*i=j0 // 0<=i<=7<br />x3=2*i=j3*i // 0<=i<=7 <br /><br /><br />fft(x1,8)的结果<br />56.0000<br />-8.0000+19.3137j<br />-8.0000+8.0000j<br />-8.0000+3.3137j<br />-8.0000<br />-8.0000-3.3137j<br />-8.0000-8.0000j<br />-8.0000-19.3137j<br />fft(x2,8)的结果<br />84.0000<br />-12.0000+28.9706j<br />-12.0000+12.0000j<br />-12.0000+4.9706j<br />-12.0000<br />-12.0000-4.9706j<br />-12.0000-12.0000j<br />-12.0000-28.9706j<br /><br />fft(x3,8)的结果经分离后的两个序列<br />56.0000<br />-8.0000+19.3137j<br />-8.0000+8.0000j<br />-8.0000+3.3137j<br />-8.0000<br />-8.0000-3.3137j<br />-8.0000-8.0000j<br />-8.0000-19.3137j<br />//******************<br />84.0000<br />-12.0000-28.9706j<br />-12.0000-12.0000j<br />-12.0000-4.9706j<br />-12.0000<br />-12.0000+4.9706j<br />-12.0000+12.0000j<br />-12.0000+28.9706j<br /><br />仿真用的程序: <br />声明:FFT主程序是在网上找的经修改后的程序,作者的相关信息均在哈<br /><br />/*时间抽选基2FFT及IFFT算法C语言实现*/<br />/*Author :Junyi Sun*/<br />/*Copyright 2004-2005*/<br />/*Mail:ccnusjy@yahoo.com.cn*/<br />#include <iostream.h><br />#include <stdio.h><br />#include <math.h><br />#include <stdlib.h><br /><br />/*定义复数类型*/<br />typedef struct<br />{<br /> float real;<br /> float img;<br />}complex;<br /><br /> void fft(complex *dt,int n); /*快速傅里叶变换*/<br />// void ifft(); /*快速傅里叶反变换*/<br /> void initW(complex *W,int n); /*初始化变换核*/<br /><br /> void depart(complex *x1,complex *x2,int n);<br /><br /> void add(complex ,complex ,complex *); /*复数加法*/<br /> void mul(complex ,complex ,complex *); /*复数乘法*/<br /> void sub(complex ,complex ,complex *); /*复数减法*/<br /> void output(complex *dt,int n); /*输出结果*/<br /><br />double PI; /*圆周率*/<br />double PI2;// 2*PI<br />int main()<br />{<br /> int i,method;<br /> int size_x=0; /*输入序列的大小，在本程序中仅限2的次幂*/<br /> double a;<br /> complex *x,*y;<br /> <br /><br /> system("cls");<br /> //PI=atan(1)*4;<br /> PI=3.141592653589793;<br /> PI2=6.283185307179586;<br /><br /> cout<<"print N=
";<br /> cin>>size_x;<br /> a=PI2/size_x;<br /><br /><br /> x=new complex;<br /> y=new complex;<br /> if(x==NULL || y==NULL)<br /> exit(0);<br /><br /> cout<<"x.real=2*i x.img=0
";<br /> for(i=0;i<size_x;i++)<br /> {<br /> x.real=2*i;<br /> x.img=0;<br /> }<br /> fft(x,size_x);<br /> output(x,size_x);<br /> <br /> cout<<"any key NO. ";<br /> cin>>i; <br /><br /> cout<<"x.real=3*i x.img=0
";<br /> for(i=0;i<size_x;i++)<br /> {<br /> x.real=3*i;<br /> x.img=0; <br /> }<br /> fft(x,size_x);<br /> output(x,size_x);<br /><br /> cout<<"any key NO. ";<br /> cin>>i; <br /><br /> cout<<"2*i+j3*i
";<br /> for(i=0;i<size_x;i++)<br /> {<br /> x.real=2*i;<br /> x.img=3*i;<br /> }<br /> fft(x,size_x); <br /> depart(x,y,size_x); //数据分离<br /> output(x,size_x);<br /> output(y,size_x);<br /><br /><br /> delete[] y;<br /> delete[] x;<br /><br /> return 0;<br />}<br /><br />/*快速傅里叶变换*/<br />void fft(complex *dt,int n)<br />{<br /> int i=0,j=0,k=0,r=0;<br /> complex up,down,product,*W;<br /><br /> double t;<br /><br />/*变址计算，将x(n)码位倒置*/<br />//***********************<br /><br /> for(i=0;i<n;i++)<br /> {<br /> k=i;j=0;<br /> t=(log(n)/log(2));<br /> while( (t--)>0 )<br /> {<br /> j=j<<1;<br /> j|=(k & 1);<br /> k=k>>1;<br /> }<br /> if(i<j)<br /> {<br /> product=dt;<br /> dt=dt[j];<br /> dt[j]=product;<br /> }<br /> }<br /><br />//*********************<br /> W=new complex[n];<br /> if(W==NULL)<br /> exit(0);<br /><br /> initW(W,n);<br />//**********************************<br /> i=j=k=0;<br /> t=(log(n)/log(2));<br /><br /> for(i=0;i<t ;i++)<br /> { //*一级蝶形运算 <br /> r=1<<i;//取 r<br /> for(j=0;j<n;j+= 2*r )<br /> { //*一组蝶形运算<br /> for(k=0;k<r;k++)<br /> { //*一个蝶形运算<br /> mul(dt[j+k+r],W[n*k/2/r],&product);<br /> add(dt[j+k],product,&up);<br /> sub(dt[j+k],product,&down);<br /> dt[j+k]=up;<br /> dt[j+k+r]=down;<br /> }<br /> }<br /> }<br /> delete[] W;<br />}<br /><br /><br />void depart(complex *x1,complex *x2,int n)<br />{<br /> float temp1_data_r,temp1_data_i,temp2_data_r,temp2_data_i;<br /> int i;<br /> complex *s_x1;<br /> <br /> <br /> s_x1=new complex[n];<br /> if(s_x1==NULL)<br /> exit(0);<br /><br /> for(i=0;i<n;i++)<br /> s_x1=x1;<br /><br /> for(i=1;i<n;i++)<br /> {<br /> temp1_data_r=(s_x1.real+s_x1[n-i].real)*0.5;<br /> temp2_data_i=(s_x1.real-s_x1[n-i].real)*0.5;<br /><br /> temp2_data_r=(s_x1.img+s_x1[n-i].img)*0.5;<br /> temp1_data_i=(s_x1.img-s_x1[n-i].img)*0.5;<br /><br /> x1.real=temp1_data_r;<br /> x1.img=temp1_data_i;<br /><br /> x2.real=temp2_data_r;<br /> x2.img=temp2_data_i;<br /> }<br /> x1[0].real=s_x1[0].real;<br /> x1[0].img=0;<br /><br /> x2[0].real=s_x1[0].img;<br /> x2[0].img=0;<br /> <br /> delete[] s_x1;<br />}<br /><br /><br /><br /><br />/*初始化变换核,cos,sin表*/<br />void initW(complex *W,int n)<br />{<br /> /*输入序列,变换核*/<br /> int i;<br /> double a;<br /><br /> a=PI2/n;<br /> for(i=0;i<n;i++)<br /> {<br /> W.real=cos(a*i);<br /> W.img=-1*sin(a*i);<br /> }<br />}<br /><br /><br />/*输出傅里叶变换的结果*/<br />void output(complex *dt,int n)<br />{<br /> int i;<br /> printf("The result are as follows
");<br /> for(i=0;i<n;i++)<br /> {<br /> printf("%.4f",dt.real);<br /><br /> if(dt.img>=0.0001)<br /> printf("+%.4fj
",dt.img);<br /> else if(fabs(dt.img)<0.0001)<br /> printf("
");<br /> else<br /> printf("%.4fj
",dt.img);<br /> }<br />}<br /><br />void add(complex a,complex b,complex *c)<br />{<br /> c->real=a.real+b.real;<br /> c->img=a.img+b.img;<br />}<br /><br />void mul(complex a,complex b,complex *c)<br />{<br /> c->real=a.real*b.real - a.img*b.img;<br /> c->img=a.real*b.img + a.img*b.real;<br />}<br /><br />void sub(complex a,complex b,complex *c)<br />{<br /> c->real=a.real-b.real;<br /> c->img=a.img-b.img;<br />}<br /><br />//可能FFT程序中还有一些BUG,请指教一下<br /><br /><br />
|
楼主
