dft.txt

数字信号处理dsp的一些c语言源代码

#include"graphics.h"
#include"math.h"
#include"stdio.h"
#include "conio.h"
#define PI 3.14159265

void figure(float x[],int m)                         /*    显示函数      */
    {  int gdriver,gmode,i; float mul=0;
          gdriver=VGA;
          gmode=VGAHI;
      initgraph(&gdriver,&gmode," ");
      setbkcolor(WHITE);
      setcolor(2);
      line(10,240,610,240);
      line(10,10,10,500);
      line(610,240,570,200);
      line(610,240,570,280);
      line(10,10,0,15);
      line(10,10,20,15);
      for(i=0;i<m;i++)
       { if(x[i]>mul)  {mul=x[i];}
       }
      mul=mul/240;
      for(i=0;i<m;i++)
      {setcolor(2);
       line(600/m*i+10,240-x[i]/mul,600/m*i+10,240);
      }
    }

void DFT(float x[],int m,float RX[],float CX[])          /*    离散傅里叶变换函数      */
   {
    int k,n;
  for(k=0;k<=m-1;k++)
     {
     RX[k]=0;CX[k]=0;
         {
          for(n=0;n<=m-1;n++)
             {
               RX[k]=RX[k]+x[n]*(cos(2*3.14*n*k/m));
               CX[k]=CX[k]+x[n]*(sin(2*3.14*n*k/m));
             }
         }
     }
    }



void IDFT(float xx[],float yx[],int m,float RX[],float CX[])/*IDFT*/
  {     int k,n;
    for(k=0;k<=m-1;k++)
     { RX[k]=0;CX[k]=0;
           {for(n=0;n<=m-1;n++)
          {RX[k]=RX[k]+xx[n]*(cos(2*PI*n*k/m))-yx[n]*(sin(2*PI*n*k/m));
           CX[k]=CX[k]+xx[n]*(sin(2*PI*n*k/m))+yx[n]*(cos(2*PI*n*k/m));
               }
        }
           RX[k]=RX[k]/m;
           CX[k]=CX[k]/m;
      }
  }


/*void IDFTR(float x[],int m,float RX[],float CX[])          反变换  没有除N
   {
    int k,n;
  for(k=0;k<=m-1;k++)
     {
     RX[k]=0;CX[k]=0;
         {
          for(n=0;n<=m-1;n++)
             {
               RX[k]=RX[k]+x[n]*(cos(2*3.14*n*k/m));
               CX[k]=CX[k]+x[n]*(-sin(2*3.14*n*k/m));
             }
         }
     }
    }

void IDFTI(float x[],int m,float RX[],float CX[])     反变换  没有除N
    {
    int k,n;
  for(k=0;k<=m-1;k++)
     {
     RX[k]=0;CX[k]=0;
         {
          for(n=0;n<=m-1;n++)
             {
               RX[k]=RX[k]+x[n]*(cos(2*3.14*n*k/m+1.57));
               CX[k]=CX[k]+x[n]*(-sin(2*3.14*n*k/m+1.57));
             }
         }
     }
    }
         */





void coin(float x[],int M, float h[], int N,float y[])     /*   卷积       */
  {
   int n,m;
   double w=0; ;
   for (n=0;n<M+N+1;n++)
   {
     w=0;
      for(m=0;m<M;m++)
      {
        if(((n-m)>=0)&&((n-m)<N))
         w+=x[m]*h[n-m];
      }
       y[n]=w;
    }
}




void main()

{
    float xn[32],hn[32],XR[32],XI[32],HR[32],HI[32],y[63],xmo[32],hmo[32],a[63]={0},b[63]={0},c[63]={0},d[63]={0},A[63]={0},B[63]={0},C[63]={0},D[63]={0};    /*,IXR[32],IXI[32],IXRR[32],IXRI[32],IXIR[32],IXII[32]; */
    float RXH[63],IXH[63],RXHR[63],RXHI[63],IXHR[63],IXHI[63],idftr[63],idfti[63],ymo[63];/*       */
    int i,j;
    for(i=0;i<32;i++)
       xn[i]=i+1;

    for(j=0;j<32;j++)
        hn[j]=100;
    figure(xn,32);   /* 显示xn     */
    getch();
    clrscr();
    figure(hn,32);         /*   显示hn       */
    getch();
    clrscr();
    DFT(xn,32,XR,XI);      /* xn 离散傅里叶变换        */
    figure(XR,32);       /*      显示变换后的实部    */
    getch();
    clrscr();
    figure(XI,32);       /*  显示变换后的虚部        */
    getch();
    clrscr();

    for(i=0;i<32;i++)
      {                              /* xn 的福利也 变换后的模植  */
       xmo[i]=sqrt(XR[i]*XR[i]+XI[i]*XI[i]);
      }
       figure(xmo,32);
       getch();
       clrscr();
    DFT(hn,32,HR,HI);      /* hn离散傅里叶变换         */
    figure(HR,32);           /*   显示变换后的实部       */
    getch();
    clrscr();
    figure(HI,32);          /*   显示变换后的虚部       */
    getch();
    clrscr();

     for(i=0;i<=32;i++)        /*显示 hn傅立叶变换后的序列*/
       {
         hmo[i]=sqrt(HR[i]*HR[i]+HI[i]*HI[i]);
       }
         figure(hmo,32);
       getch();
       clrscr();



    coin(xn,32,hn,32,y);      /*  xn   hn的卷积        */
    figure(y,63);
    getch();
    clrscr();

     for(i=0;i<=31;i++)
        {
          A[i]=xn[i];
          B[i]=hn[i];
        }

        DFT(A,63,a,b);
        DFT(B,63,c,d);






    for(i=0;i<=63;i++)                    /* 计算X乘以H后的实部与虚部 */
      {
        RXH[i]=a[i]*c[i]-b[i]*d[i];  /*  ac-bd  */
        IXH[i]=a[i]*d[i]+b[i]*c[i];   /*  ad+bc */
      }

    IDFT(RXH,IXH,63,RXHR,RXHI);      /* 对X乘H后的实部傅里叶反变换  得出的实部傅里叶反变换后的实部与虚部*/
                                     /* IDFTI(IXH,63,IXHR,IXHI); 对X乘H后的虚部傅里叶反变换 得出的虚部傅里叶反变换后的实部与虚部*/
    for(i=0;i<63;i++)
       {
          ymo[62-i

          ]=sqrt(RXHR[i]*RXHR[i]+RXHI[i]*RXHI[i]);
       }
    figure(ymo,63);
    getch();
    clrscr();





   /* for(i=0;i<63;i++)
        {
         idftr[i]=RXHR[i]+IXHR[i];
         idfti[i]=RXHI[i]+IXHI[i];
        }
    for(i=0;i<63;i++)
       {
         a[i]=a[i]*idfti[i];
       }
    figure(a,63);
    getch();
    clrscr(); */
    /*IDFTR(XR,32,IXRR,IXRI);
    figure(IXRR,32);
    getch();
    clrscr();
    figure(IXRI,32);
    getch();
    clrscr();
    IDFTR(XI,32,IXIR,IXII);
    figure(IXIR,32);
    getch();
    clrscr();
    figure(IXII,32);
    getch();
    clrscr(); */
    /*for(i=0;i<=32;i++)
      {
         xmo[i]=sqrt(pow(XR[i],2)+pow(XI[i],2));
         printf("%f\n",xmo[i]);
       }
    getch();  */

}