[分享]用matlab编的fft变换程序

francexuqu

用蝶形算法和码位倒置法编写的fft变换程序，并验证之，有兴趣的可以看看。

1. function data=myfft(datat,nn,isign)
   
2. datat=mybitrevorder(datat,nn);                %这个地方也可以用matlab自带的bitrevorder,有兴趣的还是自己编一下
   
3. for i=0:length(datat)-1
   
4.     data(2*i+1)=datat(i+1); data(2*i+2)=0;
   
5. end
   
6. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
7. n=nn.*2; mmax=2;
   
8. while n>mmax;
   
9.     istep=2.*mmax; theta=6.28318530717959/(isign*mmax); wtemp=sin(0.5*theta);
   
10.     wpr = -2.0*wtemp*wtemp; wpi=sin(theta); wr=1.0; wi=0.0;
    
11.     for m=1:2:mmax
    
12.         for i=m:istep:n
    
13.             j=i+mmax; tempr=wr*data(j)-wi*data(j+1); tempi=wr*data(j+1)+wi*data(j);
    
14.            data(j)=data(i)-tempr; data(j+1)=data(i+1)-tempi;
    
15.            data(i)=data(i)+tempr; data(i+1)=data(i+1)+tempi;
    
16.         end
    
17.         wtemp=wr; wr=wtemp*wpr-wi*wpi+wr; wi=wi*wpr+wtemp*wpi+wi;
    
18.     end
    
19.     mmax=istep;
    
20. end  
    
21. %%%%%%%%%%%%%码位倒置程序%%%%%%%%%%%%%%%%%%%%
    
22. function  A=mybitrevorder(A,N)
    
23. j=N./2;
    
24. for i=2:1:N-1
    
25.     if i<j
    
26.         temp=A(i); A(i)=A(j+1); A(j+1)=temp;
    
27.     end
    
28.     k=N/2;
    
29.     while j>=k
    
30.         j=j-k; k=k/2;
    
31.     end
    
32.     j=j+k;
    
33. end
    
34. 
    
35. %%%%%%%%%%%%%进行验证(matlab自带的fft，自己编写的fft，和函数的解析fft三者的比较%%%%%%%
    
36. clear all;
    
37. isign=1; nn=256; t=linspace(0,3,nn); f=2.*exp(-3.*t);
    
38. data=myfft(f,nn,isign);
    
39. Ts=t(2)-t(1); Ws=2.*pi./Ts; W=Ws.*(0:nn./2)./nn;
    
40. for i=0:nn-1
    
41.     datar(i+1)=data(2.*i+1).*Ts; datai(i+1)=data(2.*i+2).*Ts;
    
42. end
    
43. for i=1:nn./2+1
    
44.     co(i)=datar(i); si(i)=datai(i); power(i)=(co(i).^2+si(i).^2).^0.5;
    
45. end
    
46. F=fft(f); Fp=F(1:nn./2+1)*Ts; Fpr=real(Fp); Fpi=imag(Fp);
    
47. W=Ws.*(0:nn./2)./nn; A=abs(Fp); Fa=2./(3+j.*W); B=abs(Fa);
    
48. plot(W,B,'o',W,A,'+r',W,power,'g');

复制代码

%%%%%%结果见附件%%%%%%%%%%%

花如月

是自己写的，那可真是高手了。我都写不出来、、、、

sigma665

这里卧虎藏龙

francexuqu

其实fft的蝶形算法都已经成行好久了
C++程序都有现成的；
只要搞懂蝶形算法的来龙去脉，和码位倒置的编写就很简单了；
但是因该注意的是C++中的脚标是从0开始；
而matlab是从1开始的；
ps:本人也是刚学习matlab，希望能在这里和大家共同进步；

hyl2323

我用C++写了蝶形算法和另外一种算法，这是matlab版，能发不？

sigma665

发吧
哈哈

hyl2323

4#举重若轻，佩服佩服！可能我智商平平，我看懂原理后，也是费了好多脑力才编写出来的。晚上来发！

我以主题帖发了我用C++写的FFT，两种算法，大伙顶一顶哈，第一次发原创！

[ 本帖最后由 ChaChing 于 2010-6-25 21:16 编辑 ]

无水1324

支持原创的东西，:lol

hyl2323

无水是湖南长沙的？老乡哦，见过，握握手！

gagmeng

支持原创的东西，

cwt1128

真是太强了，很有用的东西，谢！

鹿鸣

高手如云。

ztx106

我想请问一下，如果对功率谱进行IFFT变换，再对IFFT结果进行FFT变换，那么这两种情况下的功率谱曲线是完全重合的吗？还是FFT变换结果的曲线应该有波动？

yuejunefong

看不明白 还是菜鸟一个

歪猴