与大家共享一些资料——DSP应用程序

lingyufeng

/*时间抽选基2FFT及IFFT算法C语言实现*/
/*Author :Junyi Sun*/
/*Copyright 2004-2005*/
/*Mail:ccnusjy@yahoo.com.cn*/
#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(&quotlease input the size of x:\n");
     scanf("%d",&size_x);
     printf(&quotlease 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--){

lingyufeng

很实用的工程资料 与大家共享

XZY8011

谢谢与大家共享

pedro888

thanks for your information.....
thanks....

dockiller

it is very helpful for me
thank you

icestorm

太好了，谢谢楼主！！！

tangdk

haodongxi

fengyun01

很实用的工程资料