大家好,最近我遇到了有关FFT算法的问题:我试图参考教科书在C语言中对该算法进行编程,但它不能作为matlab的fft函数使用. 首先,我在"Complex.h"中定义了用于复杂操作的数据类型.
Hello,everyone,recently I had met a problem about FFT algorithm:I tried to program the algorithm in C with reference to the textbook,but it didn''t work as the matlab'' fft function. Firstly,I defined a data type for complex operation in "Complex.h"
typedef struct { double x_real; double x_imag; }Complex; // the multiplication for complex x1 and x2 Complex CMul(Complex x1,Complex x2) { Complex result; result.x_real=x1.x_real*x2.x_real-x1.x_imag*x2.x_imag; result.x_imag=x1.x_real*x2.x_imag+x1.x_imag*x2.x_real; return result; } //the addition for the complex x1 and x2 Complex CAdd(Complex x1,Complex x2) { Complex result; result.x_real=x1.x_real+x2.x_real; result.x_imag=x1.x_imag+x2.x_imag; return result; } //the subtract for the complex x1 and x2 Complex CSub(Complex x1,Complex x2) { Complex result; result.x_real=x1.x_real-x2.x_real; result.x_imag=x1.x_imag-x2.x_imag; return result; }
然后我开始在"fft.cpp"中对FFT进行编程
Then I began to program the FFT in"fft.cpp"
#include <stdio.h> #include <math.h> #include "Complex.h" #define pi 3.1415926 //bit-reverse operation for the input void RSort(Complex x[],int length ) { int LH,J,N1,I,K; Complex T; LH=length/2; J=LH; N1=length-2; for (I=1;I<=N1;I++) { if (I<J) { T=x[I]; x[I]=x[J]; x[J]=T; } K=LH; while(J>=K) { J=J-K; K=K/2; } J=J+K; } } /* void FFT(Complex x[],int N ,int M)---compute the FFT of the origin signal x[]--the input of origin and output of its FFT N--the length of x[]; M--the power of N; */ void FFT(Complex x[],int N ,int M) { /* L--the level of the Butterfly operation B--the distance between two signal in every Butterfly operation wn[]---the rotation factor p--the exponent of rotation factor */ int L,B,J,p,k,i; Complex wn[4]; //to get the rotation factor for (i=0;i<N/2;i++) { wn[i].x_real=cos(2*pi*i/N); wn[i].x_imag=-sin(2*pi*i/N); } //bit-reversal for the input signal RSort(x,N); //the main part of the FFT algorithm for (L=1;L<=M;L++) { B=int(pow(2,L-1)); for (J=0;J<=B-1;J++) { p=int(J*pow(2,M-L)); for (k=J;k<=N-1;) { x[k]=CAdd(x[k],CMul(wn[p],x[k+B])); x[k+B]=CSub(x[k],CMul(wn[p],x[k+B])); k=k+int(pow(2,L)); } } } } void main() { Complex *x; int i; for (i=0;i<=15;i++) { x[i].x_real=i; x[i].x_imag=0; } printf("the input signal:"); for (i=0;i<=15;i++) { printf("%lf ",x[i].x_real); } printf("\n"); FFT(x,16,4); printf("the result after FFT:"); for (i=0;i<=15;i++) { printf("%lf i%lf ",x[i].x_real,x[i].x_imag); } printf("\n"); }
为方便起见,我将输入信号定义为信号x [] = [0,1,2,3,4,5,6,7].执行程序后的结果输出为:
For convenience,I defined the input signal to be a signal x[]=[0,1,2,3,4,5,6,7].The output of the result after executing program is :
报价:输入信号:0.000000 1.000000 2.000000 3.000000 4.000000 5.000000 6. 000000 7.000000 FFT后的结果:28.000000 i0.000000 -1.414213 i-4.828427 4.000000 i-6.00 0000 -0.707107 i-0.707107 12.000000 i0.000000 0.000000 i-2.000000 4.000000 i 0.000000 0.000000 i0.000000 按任意键继续
the input signal:0.000000 1.000000 2.000000 3.000000 4.000000 5.000000 6. 000000 7.000000 the result after FFT:28.000000 i0.000000 -1.414213 i-4.828427 4.000000 i-6.00 0000 -0.707107 i-0.707107 12.000000 i0.000000 0.000000 i-2.000000 4.000000 i 0.000000 0.000000 i0.000000 Press any key to continue
而且在matlab中使用FFT函数的结果是:
And the result using FFT function in matlab is:
报价:>> x = [0 1 2 3 4 5 6 7]; >> y = fft(x) y = 28.0000 -4.0000 + 9.6569i -4.0000 + 4.0000i -4.0000 + 1.6569i -4.0000 -4.0000-1.6569i -4.0000-4.0000i -4.0000-9.6569i
>> x=[0 1 2 3 4 5 6 7]; >> y=fft(x) y = 28.0000 -4.0000 + 9.6569i -4.0000 + 4.0000i -4.0000 + 1.6569i -4.0000 -4.0000 - 1.6569i -4.0000 - 4.0000i -4.0000 - 9.6569i
我非常期待您的意见,非常感谢!^ _ ^^ _ ^ ...
I''m looking forward to your opinions,Thanks very much!^_^^_^...
推荐答案请注意:您的问题并不是真的要求专家意见或类似的东西;这确实是为您做大量工作的要求:测试,调试,检测错误并尝试修复它们.也许有人会为您找到一些东西,但是对于像我这样的许多人来说,从头开始实现该算法会容易得多. 所有这些实际上都超出了CodeProject快速问答的格式.答案. 所以我不知道还有什么帮助... 至少使用C ++复杂库,而不是基本的实现. FFT有许多可以使用的实现方式.阅读一些文章,例如: 如何实现FFT算法 [ ^ ], www.librow/articles/article-10 [ ^ ]. 您可以找到更多.
作为逻辑错误的错误发生在FFT()函数中: The errors which are logical errors occur in the function of FFT(): x[k]=CAdd(x[k],CMul(wn[p],x[k+B])); x[k+B]=CSub(x[k],CMul(wn[p],x[k+B]));
这些代码应修改为(在定义Complex的这些变量T1和T2之前):
These codes should be modified to(before these variable T1 and T2 of Complex should be defined):
T2=CMul(wn[p],x[k+B]); T1=x[k]; x[k]=CAdd(T1,T2); x[k+B]=CSub(T1,T2);然后,它工作正常!对此我感到很兴奋.现在,我首先要感谢SAKryukov的帮助!:-) :-)我已经开始对自己充满信心了!再次感谢. 该程序无法有效运行,仍然需要改进!让我们继续吧!
And then,it works correctly!I felt excited about that.Now I wanted to give thanks to SAKryukov for his help firstly!:-):-)I had began to gain strong confidence in myself! Thanks again. This program didn''t work efficiently and it still needed to be improved well!Let''s go for it!
更多推荐
我的FFT算法有什么问题?
发布评论