[讨论]关于fft算法问题

simon21

傅立叶变换的程序贴出如下:

1. 
   
2. 
   
3. #include <stdio.h>
   
4. #include <math.h>
   
5. #define N 64
   
6. #define m 6 //2^m=N
   
7. /*
   
8. float twiddle[N/2]={1.0,0.707,0.0,-0.707,};
   
9. float x_r[N]={1,1,1,1,0,0,0,0,};
   
10. float x_i[N]; //N=8
    
11. */
    
12. float twiddle[N/2]={1,0.9951,0.9808,0.9570,0.9239,0.8820,0.8317,0.7733,
    
13. 0.7075,0.6349,0.5561,0.4721,0.3835,0.2912,0.1961,0.0991,
    
14. 0.0000,-0.0991,-0.1961,-0.2912,-0.3835,-0.4721,-0.5561,-0.6349,
    
15. -0.7075,-0.7733,0.8317,-0.8820,-0.9239,-0.9570,-0.9808,-0.9951,}; //N=64
    
16. float x_r[N]={1,1,1,1,1,1,1,1,
    
17. 1,1,1,1,1,1,1,1,
    
18. 1,1,1,1,1,1,1,1,
    
19. 1,1,1,1,1,1,1,1,
    
20. 0,0,0,0,0,0,0,0,
    
21. 0,0,0,0,0,0,0,0,
    
22. 0,0,0,0,0,0,0,0,
    
23. 0,0,0,0,0,0,0,0,};
    
24. float x_i[N];
    
25. 
    
26. void fft_init()
    
27. {
    
28. int i;
    
29. for(i=0;i<N;i++)
    
30. x_i[i]=0.0;
    
31. }
    
32. void bitrev()
    
33. {
    
34. int p=1,q,i;
    
35. int bit_rev[N];
    
36. float xx_r[N];
    
37. bit_rev[0]=0;
    
38. while(p<N)
    
39. {
    
40. for(q=0;q<p;q++)
    
41. {
    
42. bit_rev[q]=bit_rev[q]*2;
    
43. bit_rev[q+p]=bit_rev[q]+1;
    
44. }
    
45. p=p*2;
    
46. }
    
47. for(i=0;i<N;i++)
    
48. {
    
49. xx_r[i]=x_r[i];
    
50. }
    
51. for(i=0;i<N;i++)
    
52. {
    
53. x_r[i]=xx_r[bit_rev[i]];
    
54. }
    
55. }
    
56. 
    
57. void display()
    
58. {
    
59. int i;
    
60. for(i=0;i<N;i++)
    
61. printf("%f\t%f\n",x_r[i],x_i[i]);
    
62. }
    
63. 
    
64. void fft1()
    
65. {
    
66. int L,i,b,j,p,k,tx1,tx2;
    
67. float TI,TR,temp;
    
68. float tw1,tw2;
    
69. 
    
70. for(L=1;L<=m;L++)
    
71. { /* for(1) */
    
72. b=1; i=L-1;
    
73. while(i>0)
    
74. {b=b*2;i--;} /* b= 2^(L-1) */
    
75. for(j=0;j<=b-1;j++) /* for (2) */
    
76. { p=1; i=m-L;
    
77. while(i>0) /* p=pow(2,7-L)*j; */
    
78. {p=p*2; i--;}
    
79. p=p*j;
    
80. tx1=p%(N);
    
81. tx2=tx1+(3*N)/4;
    
82. tx2=tx2%(N);
    
83. if (tx1>=(N/2))
    
84. tw1=-twiddle[tx1-(N/2)];
    
85. else
    
86. tw1=twiddle[tx1];
    
87. if (tx2>=(N/2))
    
88. tw2=-twiddle[tx2-(N/2)];
    
89. else
    
90. tw2=twiddle[tx2];
    
91. for(k=j;k<N;k=k+2*b) /* for (3) */
    
92. {TR=x_r[k]; TI=x_i[k]; temp=x_r[k+b];
    
93. x_r[k]=x_r[k]+x_r[k+b]*tw1+x_i[k+b]*tw2;
    
94. x_i[k]=x_i[k]-x_r[k+b]*tw2+x_i[k+b]*tw1;
    
95. x_r[k+b]=TR-x_r[k+b]*tw1-x_i[k+b]*tw2;
    
96. x_i[k+b]=TI+temp*tw2-x_i[k+b]*tw1;
    
97. }
    
98. }
    
99. }
    
100. 
     
101. }
     
102. 
     
103. void dft()
     
104. {
     
105. int i,n,k,tx1,tx2;
     
106. float tw1,tw2;
     
107. float xx_r[N],xx_i[N];
     
108. //clear any data in Real and Imaginary result arrays prior to DFT
     
109. for(k=0;k<=N-1;k++)
     
110. xx_r[k]=xx_i[k]=x_i[k]=0.0;
     
111. //caculate the DFT
     
112. for(k=0;k<=(N-1);k++)
     
113. {
     
114. for(n=0;n<=(N-1);n++)
     
115. {
     
116. tx1=(n*k);
     
117. tx2=tx1+(3*N)/4;
     
118. tx1=tx1%(N);
     
119. tx2=tx2%(N);
     
120. if (tx1>=(N/2))
     
121. tw1=-twiddle[tx1-(N/2)];
     
122. else
     
123. tw1=twiddle[tx1];
     
124. if (tx2>=(N/2))
     
125. tw2=-twiddle[tx2-(N/2)];
     
126. else
     
127. tw2=twiddle[tx2];
     
128. xx_r[k]=xx_r[k]+x_r[n]*tw1;
     
129. xx_i[k]=xx_i[k]+x_r[n]*tw2;
     
130. }
     
131. xx_i[k]=-xx_i[k];
     
132. }
     
133. //display
     
134. for(i=0;i<N;i++)
     
135. printf("%f\t%f\n",xx_r[i],xx_i[i]);
     
136. 
     
137. }
     
138. 
     
139. void main()
     
140. {
     
141. /*
     
142. fft_init();
     
143. bitrev();
     
144. fft1();
     
145. display(); //FFT
     
146. */
     
147. dft(); //DFT
     
148. }
     

复制代码

simon21

　　其中分别使用fft和dft做变换

　　源程序如上: 其中函数fft1()代表fft变换,函数dft()代表dft变换.

　　当取N=8的时候,fft和dft算出来的数据是一样的.

　　fft和dft计算结果如下:

　　4.000000 0.000000

　　1.000000 -2.414000

　　0.000000 0.000000

　　1.000000 -0.414000

　　0.000000 0.000000

　　1.000000 0.414000

　　0.000000 0.000000

　　1.000000 2.414000

　　Press any key to continue

　　但是,当N=64的时候,fft和dft算出的数据x_r[ i ],x_ i[i ]就不相同:

　　fft算出的结果如下:

　　32.000000 0.000000

　　0.968460 -20.365383

　　0.000000 0.000000

　　0.993620 -6.749197

　　0.000000 0.000000

　　3.875856 -1.388196

　　0.000000 0.000000

　　... ... ... ... ;这里只取i=0~6的x_r[ i ],x_ i[i ]显示

　　dft算出的结果如下:

　　32.000000 0.000000

　　2.663399 -18.705200

　　-0.000000 0.000001

　　2.663400 -8.407199

　　0.000000 0.000000

　　2.663400 -2.332000

　　0.000000 0.000000

　　... ... ... ... ;这里只取i=0~6的x_r[ i ],x_ i[i ]显示

　　书上说,dft和fft的计算结果应该是完全一样的.这是怎么会事?是不是程序中有问题,谢谢!

Anonymous

　　fft的源程序不是现成的吗?看看这个网站，上面有世界上最优秀的fft算法，而且是免费的，matlab的fft就是用的他的算法。

　　http://www.fftw.org/

simon21

在某些场合也不完全通用吧，需要自己更具情况写

Anonymous

　　以下是引用simon21在2005-9-23 9:00:07的发言：

　　在某些场合也不完全通用吧，需要自己更具情况写

　　有啥不能通用的地方说来听听，长长见识啦，偶没有自己编过，

　　//一向都觉得只是调用个子程序罢了

　　[em06]

　　

[此贴子已经被作者于2005-9-23 15:45:38编辑过]

AaronSpark

如果只是用在信号处理里这些程序基本够了，不过我见过一篇文章将FFT和其他算法结合起来做数值分析的，用通用程序就不行了

TNC

个人感觉自己写一遍对FFT的理解程度能提高一个层次<BR>要用FFT的还是要自己写一下，不管最后用不用

陆雯

我也正在编程中.

zhyuer

好主意！