dsp.h

这是有关数字信号处理的课程的源代码

/****************************************
*filename:dsp.h
*function: this file includes two functions of fft and ifft
*****************************************/
/* define struct COMPLEX */
   typedef struct{
      float real;
	  float imag;
	  }COMPLEX;
extern void fft(COMPLEX *x,int m);
extern void ifft(COMPLEX *x,int m);
/* first function :fft,in space radix 2 decimation in frequency fft*/

void fft(COMPLEX *x,int m)
{
  static COMPLEX *w;
  static int mstore=0;
  static int n=1;


  COMPLEX u,temp,tm;
  COMPLEX *xi,*xip,*xj,*wptr;

  int i,j,k,l,le,windex;
  double arg,w_real,w_imag,wrecur_real,wrecur_imag,wtemp_real;
  if(m!=mstore)
  {
   /*  free previously allocate storage and set new m*/
       if(mstore!=0) free(w);
	   mstore=m;
	   if(m==0) return;      /* if m=0 then done */
   /*  n=2 **m=fft length*/
      n=1<<m;
	   le=n/2;
/* allocate the storage for w*/
  w=(COMPLEX *)calloc(le-1,sizeof(COMPLEX));




if(!w)
{
printf("\nUnable to allocate complex w array\n");
exit(1);

}
/*calculate the w values recursively*/
arg=4.0*atan(1.0)/le;
wrecur_real=w_real=cos(arg);
wrecur_imag=w_imag=-sin(arg);
xj=w;
for(j=1;j<le;j++)
{
 xj->real=(float)wrecur_real;
 xj->imag=(float)wrecur_imag;
 xj++;
 wtemp_real=wrecur_real *w_real-wrecur_imag *w_imag;
 wrecur_imag=wrecur_real *w_imag+wrecur_imag *w_real;
 wrecur_real=wtemp_real;
}
}
/*start fft*/
le=n;
windex=1;
for(l=0;l<m;l++)
{
 le=le/2;
 /*first iteration with no multiplies*/
   for(i=0;i<n;i=i+2*le)
     {
	   xi=x+i;
	   xip=xi+le;
	   temp.real=(xi->real+xip->real);
	   temp.imag=(xi->imag+xip->imag);
	   xip->real=(xi->real-xip->real);
	   xip->imag=(xi->imag-xip->imag);
	   *xi=temp;

	 }
/*remaining iterations use stored */
  wptr=w+windex-1;
  for(j=1;j<le;j++)
  {
  u=*wptr;
   for(i=j;i<n;i=i+2*le)
  {




xi=x+i;
xip=xi+le;
temp.real=(xi->real+xip->real);
temp.imag=(xi->imag+xip->imag);
tm.real=xi->real-xip->real;
tm.imag=xi->imag-xip->imag;
xip->real=(tm.real*u.real-tm.imag*u.imag);
xip->imag=(tm.real*u.imag+tm.imag*u.real);
*xi=temp;
}
wptr=wptr+windex;
}
windex=2*windex;
}

for(i=0;i<n;++i)
{
	j=0;
   for(k=0;k<m;++k)
   j=(j<<1)|(1&(i>>k));
   if(i<j)
   { xi=x+i;
     xj=x+j;
	 temp=*xj;
	 *xj=*xi;
	 *xi=temp;
   }
}
}

void ifft(COMPLEX*x,int m)
{
   static COMPLEX *w;
   static int mstore=0;
   static int n=1;

   COMPLEX u,temp,tm;
   COMPLEX *xi,*xip,*xj,*wptr;

   int i,j,k,l,le,windex;
   float scale;
   double arg,w_real,w_imag,wrecur_real,wrecur_imag,wtemp_real;






if(m!=mstore)	
{
if(mstore!=0)free(w);
mstore=m;
if(m==0)return;
n=1<<m;
le=n/2;
w=(COMPLEX *)calloc(le-1,sizeof(COMPLEX));
if(!w)
{
printf("\nUnable to alllocate complex w array\n");
exit(1);
}
arg=4.0*atan(1.0)/le;
wrecur_real=w_real=cos(arg);
wrecur_imag=w_imag=sin(arg);
xj=w;
for(j=1;j<le;j++)
{
xj->real=(float)wrecur_real;
xj->imag=(float)wrecur_imag;
xj++;
wtemp_real=wrecur_real*w_real-wrecur_imag*w_imag;
wrecur_imag=wrecur_real*w_imag+wrecur_imag*w_real;
wrecur_real=wtemp_real;
}
}
le=n;
windex=1;
for(l=0;l<m;l++)
{
le=le/2;
for(i=0;i<n;i=i+2*le)
{
xi=x+i;
xip=xi+le;
temp.real=(xi->real+xip->real);
temp.imag=(xi->imag+xip->imag);
xip->real=(xi->real-xip->real);
xip->imag=(xi->imag-xip->imag);
*xi=temp;

}

/* remaining iterations use stored w */
  wptr=w+windex-1;
  for(j=1;j<le;j++)
   {
     u=*wptr;
     for(i=j;i<n;i=i+2*le)
       {
         xi=x+i;
         xip=xi+le;
         temp.real=(xi->real+xip->real);
         temp.imag=(xi->imag+xip->imag);
         tm.real=xi->real-xip->real;
         tm.imag=xi->imag-xip->imag;
         xip->real=(tm.real*u.real-tm.imag*u.imag);
         xip->imag=(tm.real*u.imag+tm.imag*u.real);
         *xi=temp;
        }
     wptr=wptr+windex;
   }
   windex=2*windex;
}


   for(i=0;i<n;++i)
    { 
      j=0;
      for(k=0;k<m;++k)
          j=(j<<1)|(1&(i>>k));
          if(i<j)
          {
             xi=x+i;
             xj=x+j;
             temp=*xj;
             *xj=*xi;
             *xi=temp;
          }
     }

     scale=(float)(1.0/n);
     for(i=0;i<n;i++)
         {
           x[i].real=scale*x[i].real;
           x[i].imag=scale*x[i].imag;
         }  
}