cor01.c

这是DSP编写的自相关的程序

#include "math.h"
#define PI  3.1415927
#define Length  128           /*input array x,y length*/
#define Lengthcor  2*Length-1   /*ouput array r length*/

int i,k,j,mode;
double sum_real, t_real, temp_real;
double x_real[Length],y_real[Length];
double r_real[Lengthcor];   

main()
{    
   k=0;
   sum_real=0;
   t_real=0;
   temp_real=0;
   
/* 
	when mode=1, result is biased estimate;
   	when mode=0, result is unbiased estimate
*/                               
   
   mode=1;
   
   for(i=0; i<Length-1; i++)   /*Initialize*/
   {
      x_real[i]=0;
      y_real[i]=0;
      r_real[i]=0;
   }

   for(i=0; i<=Length-1; i++)  /*Input x,y*/
   {
      x_real[i]=sin(2*PI*i/(Length-1));
      y_real[i]=sin(2*PI*i/(Length-1));
   }

   for(k=0; k<=Length-1; k++)  /* Caculation */
   {
   	  sum_real=0;
   	  
   	  for(j=0; j<=Length-1-k; j++)
      {
         t_real=x_real[j+k]*y_real[j];
         sum_real=sum_real+t_real;
      }

      if(mode==0)
      {
         r_real[Length-1-k]=sum_real/(Length-k);
      }

      else
      {
         r_real[Length-1-k]=sum_real/Length;
      }
   }
   for(k=0; k<=Length-1; k++)  
   {
      sum_real=0;

      for(j=0; j<=Length-1-k; j++)
      {
         t_real=x_real[j]*y_real[j+k];
         sum_real=sum_real+t_real;
      }

      if(mode==0)
      {
         r_real[Length-1+k]=sum_real/(Length-k);
      }

      else
      {
         r_real[Length-1+k]=sum_real/Length;
      }
   }
   for (;;)
   {
     i=0;
   }
}