使用matlab编写fft函数

Spire

闲来无事，自己写了一个fft函数的程序，写完之后与matlab自带的fft进行比较，结果不正确，不知道我的程序里面出现了什么问题，请高手指点一二。程序在附件里。

Spire

居然这么久都没有人来看这个帖子。
这个问题已经解决了。

ChaChing

怎没人来看这个帖子? 不是破百了!:@)
而且LZ的副档应该直接贴, 不然我的系统开启会乱码
还有可以的话建议LZ与大家分享你的成果, 做个结束!

[ 本帖最后由 ChaChing 于 2009-10-12 23:01 编辑 ]

Spire

贴出程序，希望对大家有帮助。不过人们都用自带的fft函数，谁会用这个呢，哈

1. 
   
2. clear; clc;
   
3. 
   
4. for kk=1:1024; %%生成一个含有1024个数的数组data（正式使用的时候load文件）
   
5.    data(kk)=kk;
   
6. end
   
7. 
   
8. 
   
9. N=1024; %%数组总元素个数
   
10. n=10; %%2的幂指数，也就是2^10=1024，n也等于快速傅里叶变换时的总行数
    
11. j=1; %%初始化j，方便在变换数组元素次序的时候使用
    
12. 
    
13. for ii=1:N %%这个for循环的结果使原来数组次序改变，如数组1 2 3 ……14 15 16经过此循环以后变为
    
14. %%以下次序：1 9 5 13 3 11 7 15 2 10 6 14 4 12 8 16
    
15. if j>ii
    
16.    temp=data(j); data(j)=data(ii); data(ii)=temp;
    
17. end
    
18. m=N/2;
    
19. while m>=2 && j>m-1
    
20.    j=j-m; m=m/2;
    
21. end
    
22. j=j+m;
    
23. end
    
24. 
    
25. a=data; %%将变换次序以后的数组赋值给a
    
26. 
    
27. %%以下是蝴蝶算法
    
28. for r=1:n %%最外层循环是蝴蝶算法的总行数
    
29. M=2^(n-r); %%M为第r行的组数（比如第一行有2^(n-1)组，第二行有2^(n-2)组……最后一行只有一组）
    
30. for m=1:M %%m为当前计算所在的行数
    
31. %%以下两个for循环处于嵌套循环的最内层，在同一组内，前一半的计算公式和后一半的计算公式不一样
    
32. for j=1:(2^(r-1))
    
33. W=cos(pi*(j-1)/2^(r-1))-i*sin(pi*(j-1)/2^(r-1));
    
34. b(j+(m-1)*(2^(r)))=a(j+(m-1)*(2^(r)))+W*a(j+(2^(r-1))+(m-1)*(2^(r)));
    
35. end
    
36. for j=(2^(r-1)+1):(2^r)
    
37. W=-cos(pi*(j-1-2^(r-1))/2^(r-1))+i*sin(pi*(j-1-2^(r-1))/2^(r-1));
    
38. b(j+(m-1)*(2^(r)))=a(j-(2^(r-1))+(m-1)*(2^(r)))+W*a(j+(m-1)*(2^(r)));
    
39. end
    
40. end
    
41. %%用新的数组b覆盖原来的数组a，用于下一行的计算
    
42. for j=1:N, a(j)=b(j); end
    
43. end
    

复制代码

[ 本帖最后由 ChaChing 于 2009-10-13 11:57 编辑 ]

xushuai0794

我挺欣赏楼住的，我喜欢你这种追根就地的性格，这才是真正的学习:@D :@D

ChaChing

其实代码可以优化下, 如
for kk=1:1024, data(kk)=kk; end
可改为data=[1:1024]; ...

Spire

谢谢楼上二位！