fft算法.txt

我收集的fft算法

http://www.vchelp.net/cndevforum/subject_view.asp?subject_id=105642&forum_id=47




正在学数字信号处理,感觉上学期信号与系统学得不扎实，因为当时只是死记公式，这学期数信老师提倡动手实践，觉得自己在编程中对公式理解得更加深刻了。
以下是我写的FFT，欢迎指教。


/*时间抽选基2FFT及IFFT算法C语言实现*/
/*Author :Junyi Sun*/
/*Copyright 2004-2005*/
/*Mail:[e-mail]*/
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define N 1000
/*定义复数类型*/
typedef struct{
     double real;
     double img;
}complex;

complex x[N], *W; /*输入序列,变换核*/
int size_x=0;     /*输入序列的大小，在本程序中仅限2的次幂*/
double PI;        /*圆周率*/

int main(){
     int i,method;
     void fft();    /*快速傅里叶变换*/
     void ifft();
     void initW();  /*初始化变换核*/
     void change(); /*变址*/
     void add(complex a,complex b,complex *c); /*复数加法*/
     void mul(complex a,complex b,complex *c); /*复数乘法*/
     void sub(complex a,complex b,complex *c); /*复数减法*/
     void divi(complex a,complex b,complex *c);/*复数除法*/
     void output();                            /*输出结果*/
     system("cls");
     PI=atan(1)*4;
     printf("Please input the size of x:\n");
     scanf("%d",&size_x);
     printf("Please input the data in x[N]:\n");
     for(i=0;i<size_x;i++)
           scanf("%lf%lf",&x.real,&x.img);
     initW();
     printf("Use FFT(0) or IFFT(1)?\n");
     scanf("%d",&method);
     if(method==0)
           fft();
     else
           ifft();
     output();
     return 0;
}

/*快速傅里叶变换*/
void fft(){
     int i=0,j=0,k=0,l=0;
     complex up,down;
     change();
     for(i=0;i< (int)( log(size_x)/log(2) );i++){  /*一级蝶形运算*/
           l=( 1<<i );
           for(j=0;j<size_x;j+= (1<<l) ){            /*一组蝶形运算*/
                 for(k=0;k<l;k++){                                /*一个蝶形运算*/
                        mul(x[j+k+l],W[size_x*k/2/l],&up);
                        add(x[j+k],up,&up);
                        mul(x[j+k+l],W[size_x*k/2/l],&down);
                        sub(x[j+k],down,&down);
                        x[j+k]=up;
                        x[j+k+l]=down;
                 }
           }
     }
}

/*快速傅里叶逆变换*/
void ifft(){
     int i=0,j=0,k=0,l=size_x;
     complex up,down;
     for(i=0;i< (int)( log(size_x)/log(2) );i++){  /*一级蝶形运算*/
           l/=2;
           for(j=0;j<size_x;j+= (1<<l) ){            /*一组蝶形运算*/
                 for(k=0;k<l;k++){                                /*一个蝶形运算*/
                       add(x[j+k],x[j+k+l],&up);
                       up.real/=2;up.img/=2;
                       sub(x[j+k],x[j+k+l],&down);
                       down.real/=2;down.img/=2;
                       divi(down,W[size_x*k/2/l],&down);
                       x[j+k]=up;
                       x[j+k+l]=down;
                 }
           }
     }
     change();
}

/*初始化变换核*/
void initW(){
     int i;
     W=(complex *)malloc(sizeof(complex) * size_x);
     for(i=0;i<size_x;i++){
           W.real=cos(2*PI/size_x*i);
           W.img=-1*sin(2*PI/size_x*i);
     }
}

/*变址计算，将x(n)码位倒置*/
void change(){
     complex temp;
     int i=0,j=0,k=0,t;

     for(i=0;i<size_x;i++){
           k=i;j=0;
           t=(unsigned) (log(size_x)/log(2));
           while(t--){
                 j=j<<1;
                 j|=(k & 1);
                 k=k>>1;
           }
           if(j>i){
                 temp=x;
                 x=x[j];
                 x[j]=temp;
           }
     }
}

/*输出傅里叶变换的结果*/
void output(){
     int i;
     printf("The result are as follows\n");
     for(i=0;i<size_x;i++){
           printf("%.4f",x.real);
           if(x.img>=0.0001)printf("+%.4fj\n",x.img);
           else if(fabs(x.img)<0.0001)printf("\n");
           else printf("%.4fj\n",x.img);
     }
}
void add(complex a,complex b,complex *c){
     c->real=a.real+b.real;
     c->img=a.img+b.img;
}

void mul(complex a,complex b,complex *c){
     c->real=a.real*b.real - a.img*b.img;
     c->img=a.real*b.img + a.img*b.real;
}
void sub(complex a,complex b,complex *c){
     c->real=a.real-b.real;
     c->img=a.img-b.img;
}
void divi(complex a,complex b,complex *c){
     c->real=( a.real*b.real+a.img*b.img )/( b.real*b.real+b.img*b.img);
     c->img=( a.img*b.real-a.real*b.img)/(b.real*b.real+b.img*b.img);
}


比如求fft([1,2j,3,4j])
运行过程：
Please input the size of x:
4
Please input the data in x[N]:
1 0
0 2
3 0
0 4
Use FFT(0) or IFFT(1)?
0
The result are as follows
4.0000+6.0000j
-4.0000
4.0000-6.0000j
0.0000










#include <stdio.h>
#include <math.h>
#define N  64
#define m  6  //2^m=N
/*
float twiddle[N/2]={1.0,0.707,0.0,-0.707,};
float x_r[N]={1,1,1,1,0,0,0,0,};
float x_i[N];            //N=8
*/
float twiddle[N/2]={1,0.9951,0.9808,0.9570,0.9239,0.8820,0.8317,0.7733,
                  0.7075,0.6349,0.5561,0.4721,0.3835,0.2912,0.1961,0.0991,
                  0.0000,-0.0991,-0.1961,-0.2912,-0.3835,-0.4721,-0.5561,-0.6349,
                  -0.7075,-0.7733,-0.8317,-0.8820,-0.9239,-0.9570,-0.9808,-0.9951,};   //N=64
float x_r[N]={1,1,1,1,1,1,1,1,
              1,1,1,1,1,1,1,1,
              1,1,1,1,1,1,1,1,
              1,1,1,1,1,1,1,1,
                 0,0,0,0,0,0,0,0,
              0,0,0,0,0,0,0,0,
              0,0,0,0,0,0,0,0,
              0,0,0,0,0,0,0,0,};
float x_i[N];

void fft_init()
{
    int i;
    for(i=0;i<N;i++)
        x_i[i]=0.0;
}
void bitrev()
{
    int p=1,q,i;
    int bit_rev[N];
    float xx_r[N];
    bit_rev[0]=0;
    while(p<N)
    {
        for(q=0;q<p;q++)
        {
            bit_rev[q]=bit_rev[q]*2;
            bit_rev[q+p]=bit_rev[q]+1;
        }
        p=p*2;
    }
    for(i=0;i<N;i++)
    {
        xx_r[i]=x_r[i];
    }
    for(i=0;i<N;i++)
    {
        x_r[i]=xx_r[bit_rev[i]];
    }
}

void display()
{
    int i;
    for(i=0;i<N;i++)
        printf("%f\t%f\n",x_r[i],x_i[i]);
}

void fft1()
{
    int L,i,b,j,p,k,tx1,tx2;
    float TI,TR,temp;
    float tw1,tw2;

    for(L=1;L<=m;L++)
    {                                 /* for(1) */
        b=1; i=L-1;
        while(i>0) 
        {b=b*2;i--;}                 /* b= 2^(L-1) */        
        for(j=0;j<=b-1;j++)             /* for (2) */
        { p=1; i=m-L;
        while(i>0)                     /* p=pow(2,7-L)*j; */
        {p=p*2; i--;}
        p=p*j;
        tx1=p%(N);
        tx2=tx1+(3*N)/4;
        tx2=tx2%(N);
        if (tx1>=(N/2))
                tw1=-twiddle[tx1-(N/2)];
            else
                tw1=twiddle[tx1];
        if (tx2>=(N/2))
                tw2=-twiddle[tx2-(N/2)];
            else
                tw2=twiddle[tx2];
        for(k=j;k<N;k=k+2*b)        /* for (3) */
        {TR=x_r[k]; TI=x_i[k]; temp=x_r[k+b];
        x_r[k]=x_r[k]+x_r[k+b]*tw1+x_i[k+b]*tw2;
        x_i[k]=x_i[k]-x_r[k+b]*tw2+x_i[k+b]*tw1;
        x_r[k+b]=TR-x_r[k+b]*tw1-x_i[k+b]*tw2;
        x_i[k+b]=TI+temp*tw2-x_i[k+b]*tw1;
        }
        }
    }

}

void dft()
{
    int i,n,k,tx1,tx2;
    float tw1,tw2;
    float xx_r[N],xx_i[N];
    //clear any data in Real and Imaginary result arrays prior to DFT
    for(k=0;k<=N-1;k++)
        xx_r[k]=xx_i[k]=x_i[k]=0.0;
    //caculate the DFT
    for(k=0;k<=(N-1);k++)
    {
        for(n=0;n<=(N-1);n++)
        {
        tx1=(n*k);
        tx2=tx1+(3*N)/4;
        tx1=tx1%(N);
        tx2=tx2%(N);
        if (tx1>=(N/2))
                tw1=-twiddle[tx1-(N/2)];
            else
                tw1=twiddle[tx1];
        if (tx2>=(N/2))
                tw2=-twiddle[tx2-(N/2)];
            else
                tw2=twiddle[tx2];
        xx_r[k]=xx_r[k]+x_r[n]*tw1;
        xx_i[k]=xx_i[k]+x_r[n]*tw2;
        }
        xx_i[k]=-xx_i[k];
    }
    //display
     for(i=0;i<N;i++)
        printf("%f\t%f\n",xx_r[i],xx_i[i]);

}

void main()
{
/*
    fft_init();
    bitrev();
    fft1();    
    display();       //FFT
*/
    dft();           //DFT       
}
[/code]
大家好,我写了一个傅立叶变换的C程序.其中分别使用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的计算结果应该是完全一样的.这是怎么会事?是不是我的程序中还有哪些问题,请傅立叶达人指点指点,谢谢.


 
不好意思,
float twiddle[N/2]={1,0.9951... ...
第四行,第三列的0.8317应该是-0.8317 呵呵,忘记加负号了.sorry
我准备这个来做谐波分析的.
程序基本上是没有问题了.希望对做傅立叶变换的朋友有帮助^_^  




以下是从网上DOWNLOAD的FFT源代码，拜托高手帮我把相应的公式及讲解写出来，最好能提供相关的文字资料，信号处理书上只是简单的提了一下，不是很明白。
/*************************************************************************
*
* 函数名称：
*   FFT()
*
* 参数:
*   complex<float> * TD    - 指向时域数组的指针
*   complex<float> * FD    - 指向频域数组的指针
*   r                        －2的幂数，即迭代次数
*
* 返回值:
*   无。
*
* 说明:
*   该函数用来实现快速付立叶变换。
*
************************************************************************/
VOID WINAPI FFT( complex<float> * TD, complex<float> * FD, int r)
{
    // 付立叶变换点数
    LONG    count;
    
    // 循环变量
    int        i,j,k;