test13.c

这里面有很多dsp应用程序源码

#include <math.h>
#include <stdio.h>

#define pi 3.1415925
int n=8;
int N=16;                   /* N=2*n */
int r=4;
int k,i,j,bfsize,p,t;
float temp1,temp2,u_real,u_imag,v_real,v_imag;
float wk_real[50],wk_imag[50];
float y_real[50],y_imag[50],z_imag[50],z_real[50];
float x_real[50],x_imag[50];
float xc[50];
float x1_real[50],x1_imag[50],w_real[50],w_imag[50];
float y1_real[50],y1_imag[50],x2_real[50],x2_imag[50];

main()
{
   for(i=0;i<n;i++)
   {
      x1_real[i]=exp(-i);
      x1_imag[i]=0;
   }

   for(i=n;i<N;i++)
   {
      x1_real[i]=0;
      x1_imag[i]=0;
   }

   for(i=0;i<N/2;i++)                  /*FFT*/
   {
      wk_real[i]=cos(-2*pi*i/N);
      wk_imag[i]=sin(-2*pi*i/N);
   }

   for(i=0;i<N;i++)
   {
      y_real[i]=x1_real[i];
      y_imag[i]=x1_imag[i];
      z_real[i]=x2_real[i];
      z_imag[i]=x2_imag[i];
   }

   for(k=0;k<r;k++)
   {
      for(j=0;j<1<<k;j++)
      {
         bfsize=1<<(r-k);
         for(i=0;i<bfsize/2;i++)
         {
            p=j*bfsize;
            t=i*(1<<k);

            u_real=y_real[i+p]-y_real[i+p+bfsize/2];
            u_imag=y_imag[i+p]-y_imag[i+p+bfsize/2];

            v_real=u_real*wk_real[t]-u_imag*wk_imag[t];
            v_imag=u_real*wk_imag[t]+u_imag*wk_real[t];

            z_real[i+p]=y_real[i+p]+y_real[i+p+bfsize/2];
            z_imag[i+p]=y_imag[i+p]+y_imag[i+p+bfsize/2];

            z_real[i+p+bfsize/2]=v_real;
            z_imag[i+p+bfsize/2]=v_imag;
         }
      }

      for(i=0;i<N;i++)
      {
         x_real[i]=y_real[i];
         x_imag[i]=y_imag[i];
	     y_real[i]=z_real[i];
         y_imag[i]=z_imag[i];
         z_real[i]=x_real[i];
         z_imag[i]=x_imag[i];
      }
   }

   for(j=0;j<N;j++)
   {
      p=0;

      for(i=0;i<r;i++)
      {
         if(j&(1<<i))
            p+=1<<(r-i-1);
      }

      x2_real[j]=y_real[p];
      x2_imag[j]=y_imag[p];                        
   }

   for(k=0;k<2*n;k++)
	{
		  w_real[k]=cos(pi*k/(2*n));
		  w_imag[k]=-sin(pi*k/(2*n));
		  y1_real[k]=x2_real[k]*w_real[k]-x2_imag[k]*w_imag[k];
		  y1_imag[k]=x2_real[k]*w_imag[k]+x2_imag[k]*w_real[k];
	}
	temp1=sqrt(1./n);
	temp2=sqrt(2./n);
	xc[0]=temp1*(y1_real[0]);
	for(k=1;k<n;k++)
	{
		   xc[k]=temp2*(y1_real[k]);
	}
}