c_fft.c

数字信号处理资料 ?中藕糯碜柿?

/*基二频域算法*/
#include "math.h"
#include "stdio.h"

struct compx 
{ double real;
  double imag;
} compx ;


struct compx EE(struct compx b1,struct compx b2)
{
struct compx b3;
b3.real=b1.real*b2.real-b1.imag*b2.imag;
b3.imag=b1.real*b2.imag+b1.imag*b2.real;
return(b3);
}

void FFT(struct compx *xin,int N)
{
int f,m,LH,nm,i,k,j,L;
double p , ps ;
int le,B,ip;
float pi;
struct compx v,w,t;
LH=N/2;  
f=N;
for(m=1;(f=f/2)!=1;m++){;}  /*2^m=N*/

{
for(L=m;L>=1;L--)    /*这里和时域的也有差别*/
{  
le=pow(2,L);
B=le/2; /*每一级碟形运算间隔的点数*/
pi=3.14159;
 for(j=0;j<=B-1;j++)
  {
   p=pow(2,m-L)*j;
   ps=2*pi/N*p;
   w.real=cos(ps);
   w.imag=-sin(ps);
   for(i=j;i<=N-1;i=i+le)
     { 
      ip=i+B;   
      t=xin[i];
      xin[i].real=xin[i].real+xin[ip].real;
      xin[i].imag=xin[i].imag+xin[ip].imag;   
      xin[ip].real=xin[ip].real-t.real;
      xin[ip].imag=xin[ip].imag-t.imag;      
      xin[ip]=EE(xin[ip],w);
     }
  }
}
} 
/*变址运算*/

nm=N-2;    
j=N/2;
for(i=1;i<=nm;i++)
{
if(i<j){t=xin[j];xin[j]=xin[i];xin[i]=t;}
k=LH;
while(j>=k){j=j-k;k=k/2;}
j=j+k;
}

} 
/******************************************************/


/*****************main programe********************/




#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#define   Num  32
float  result[Num];
struct  compx s[Num];
const float  sin_[]={0,0.19509032, 0.38268343,
0.55557023,0.70710678,0.83146961,
0.92387953,0.98078528,1,0.98078528,
0.92387953,0.83146961,0.70710678,
0.55557023,
0.38268343,0.19509032,0,
   -0.19509032, -0.38268343,
-0.55557023,-0.70710678,-0.83146961,
-0.92387953,-0.98078528,-1,-0.98078528,
-0.92387953,-0.83146961,-0.70710678,
-0.55557023,
-0.38268343,-0.19509032,
};
const float pp=3.14159;
float  hanming_table[Num];
main()
{

int i;
/*
for(i=0;i<Num;i++)
   hanming_table[i]=0.54-0.46*(2*pp*i/(Num-1));
 
*/
/*
for(i=0;i<16;i++)
{

s[i].real=sin(pp*i*2/16);


s[i].imag=0;
}
*/
 for(i=0;i<32;i++)
{

  switch(i)
  {
     case 0:s[i].real=2;break;

     case 1:s[i].real=1398;break;
     case 2:s[i].real=1294;break;

     case 3:s[i].real=1321;break;

     case 4:s[i].real=1478;break;
     case 5:s[i].real=1740;break;
case 6:s[i].real=2061;break;
case 7:s[i].real=2384;break;
   case 8:s[i].real=2667;break;
     case 9:s[i].real=2871;break;
          case 10:s[i].real=2956;break;
          case 11:s[i].real=2903;break;

   case 12:s[i].real=2729;break;
     case 13:s[i].real=2471;break;
    case 14:s[i].real=2164;break;
    case 15:s[i].real=1846;break;

    case 16:s[i].real=1565;break;

     case 17:s[i].real=1368;break;
case 18:s[i].real=1289;break;
  case 19:s[i].real=1340;break;
    case 20:s[i].real=1521;break;
    case 21:s[i].real=1799;break;
    case 22:s[i].real=2121;break;

      case 23:s[i].real=2446;break;

          case 24:s[i].real=2715;break;

          case 25:s[i].real=2896;break;
          case 26:s[i].real=2958;break;
              case 27:s[i].real=2877;break;

               case 28:s[i].real=2684;break;
               case 29:s[i].real=2416;break;
               case 30:s[i].real=2105;break;
                  case 31:s[i].real=1791;break;


  }
  s[i].real=s[i].real/1023.75-2.0+1;

s[i].imag=0;
}
/*
for(i=0;i<7;i++)
{
s[i].real=1;
s[i].imag=0;
}
for(i=8;i<15;i++)
{
s[i].real=0;
s[i].imag=0;
}
*/
/*
for(i=0;i<Num;i++)
{
s[i].real*=hanming_table[i];
s[i].imag*=hanming_table[i];
}

*/

FFT(s,Num);
/*
for(i=0;i<Num;i++)
{
s[i].real/=Num;
s[i].imag/=Num;
}
for(i=0;i<Num;i++)
{
s[i].real/=hanming_table[i];
s[i].imag/=hanming_table[i];
}
 */
for(i=0;i<16;i++)
{
printf("%.4f",s[i].real);
printf("+%.4fj\n",s[i].imag);
result[i]=sqrt(pow(s[i].real,2)+pow(s[i].imag,2));
printf("+%.4fj\n",result[i]);
     getch();
}

}