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

lingyufeng · 发表于 2006-10-25 12:26:37

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

您需要登录才可以下载或查看，没有账号？注册

x

本帖最后由 cjsb37 于 2013-4-29 09:11 编辑

/*时间抽选基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(&quot

lease input the size of x:\n");
scanf("%d",&size_x);
printf(&quot

lease 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--){