fast convolution.c

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

#include<stdio.h>
#include<stdlib.h>
#include<dos.h>
#include<math.h>
#include<conio.h>
#include"dsp.h"

void main()
{
   COMPLEX  *x,*h,*y,*xtemp,*htemp;
   int m,n,i,nx,nh;
   float  mm,mp;
   system("cls");
   printf("Input the points of x(n)......nx=");
   scanf("%d",&nx);
   printf("Input the points of h(n)......nh=");
   scanf("%d",&nh);
   /*calculate  the total points n for implementing  fast convolution*/
   n=nx+nh-1;
   mm=log10((double)n)/log10(2.0);  /*mm=log2(n)*/
   mp=mm-(int)mm;
   m=(int)mm;
   if(mp>0)     /*mp>0 means n=2 **m ,m isn't a integer*/
     {
      m=(int)mm+1;
      n=1<<m;
     }
   /*allocate complex x,h and y array*/
   x=(COMPLEX *)calloc(n,sizeof(COMPLEX));
   h=(COMPLEX *)calloc(n,sizeof(COMPLEX));
   y=(COMPLEX *)calloc(n,sizeof(COMPLEX));
   xtemp=(COMPLEX *)calloc(n,sizeof(COMPLEX));
   htemp=(COMPLEX *)calloc(n,sizeof(COMPLEX));
   
   if(!x||!h||!y)
   {
     printf("Unable to allcoate complex x,h or y array\n");
     exit(1);
   } 
   /* store x(n)=exp(-n) and h(n)=1*/
   /*this part will be changed as your need!*/
     for(i=0;i<nx;i++)
      {
        x[i].real=xtemp[i].real=(float)exp((double)(-i));
        x[i].imag=xtemp[i].imag=0.0;
        
      }
     for(i=0;i<nh;i++)
       {
       h[i].real=htemp[i].real=1.0;
       h[i].imag=xtemp[i].imag=0.0;
       
       } 
   /*add values 0 to the tail of x(n),h(n)*/
   for(i=nx;i<n;i++)x[i].real=x[i].imag=xtemp[i].real=xtemp[i].imag=0.0;
   for(i=nh;i<n;i++)h[i].real=h[i].imag=htemp[i].real=htemp[i].imag=0.0;
   /*fft for x(n),h(n)*/
     fft(x,m);
     fft(h,m);
   /*Y(K)=X(K)*H(K)*/
     for(i=0;i<n;i++)
        {
          y[i].real=x[i].real*h[i].real-x[i].imag*h[i].imag;
          y[i].imag=x[i].real*h[i].imag+x[i].imag*h[i].real;
        }
   /*IFFT for Y(K)*/
   ifft(y,m);
   /*output y(n)*/
   printf("\n");
   printf("     n     y.Real      y.Imag");
   printf("\n-------------------------------------------------\n");
   for(i=0;i<n;i++)
     {
      printf("\n %4d   %9f    %9f",i,y[i].real,y[i].imag);
     }
     printf("\n");
}