frequency_analysis.c

用FFT对连续信号和时域离散信号进行谱分析

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<conio.h>
#include<graphics.h>
void InitGraphic(void)
{int driver=DETECT,mode=0;
initgraph(&driver,&mode,"c:\\tc");
setfillstyle(1,BLACK);
bar(0,0,639,479);
return;
}
void PlotTimeSequ(float *x,int n,float scale,int space,int y1)
{int j,j1,x0,x1;
float max=0.0;
x0=(640-n*space)/2;
setcolor(WHITE);
line(x0-20,y1,640-x0+20,y1);
line(x0,y1-90,x0,y1+90);
outtextxy(x0-30,y1-80,"x(n)");
outtextxy(x0-10,y1+10,"0");
outtextxy(640-x0+5,y1-10,"n");
outtextxy(640-x0-4,y1+6,"N");
outtextxy(320-10,y1+10,"N/2");
setfillstyle(1,4);
setcolor(6);
for(j=0;j<n;j++)
    {if(x[j]>max) max=x[j];
    }
for(j=0;j<n;j++)
    {
     j1=space*j+x0;
     x1=x[j]*scale/max;
     line(j1,y1,j1,y1-x1);
     fillellipse(j1,y1-x1,2,2);
    }
 return;
}
void PlotSpectrum(float *x,int n,float scale,int space,int y1)
{
 int j,j1,x0,x1;
 float max=0.0;
 x0=(640-n*space)/2;

 setcolor(WHITE);
 line(x0-20,y1,640-x0+20,y1);
 line(x0,y1-150,x0,y1+10);
 outtextxy(x0-30,y1-140,"X(k)");
 outtextxy(x0-10,y1+10,"0");
 outtextxy(640-x0+5,y1-10,"k");
 outtextxy(640-x0-4,y1+6,"N");
 outtextxy(320-10,y1+10,"N/2");

 setfillstyle(1,4);
 setcolor(6);
 for(j=0;j<n;j++)
 {if(x[j]>max) max=x[j];
 }
 for(j=0;j<n;j++)
      {j1=space*j+x0;
        x1=x[j]*scale/max;
        line(j1,y1,j1,y1-x1);
        fillellipse(j1,y1-x1,2,2);
      }
 return;
}
void tra(float *x,float *y);
void FFT(float *xr,float *xi,int n,int inv)
{
 int i,j,a,b,k,m;
 int ep,arg,mt,s0,s1;
 float sign,pr,pi,ph;
 float *c,*s;
c=(float*)calloc(n,sizeof(float));
 if(c==NULL)exit(1);
 s=(float*)calloc(n,sizeof(float));
 if(s==NULL)exit(1);
 j=0;
 if(inv==0)
 {sign=1.0;
  for(i=0;i<n;i++)
  {xr[i]=xr[i]/n;
    xi[i]=xi[i]/n;
  }
 }
 else sign=-1.0;
 for(i=0;i<n-1;i++)
     {if(i<j)
      {tra(&xr[i],&xr[j]);
        tra(&xi[i],&xi[j]);
      }
      k=n/2;
      while(k<=j)
              {j=j-k;
            	k=k/2;
        	  }
          j=j+k;
        	 }
 ep=0;
 i=n;
 while(i!=1)
 {ep=ep+1;
  i=i/2;
 }
 ph=2*M_PI/n;

 for(i=0;i<n;++i)
      {s[i]=sign*sin(ph*i);
        c[i]=cos(ph*i);
      }
 a=2;
 b=1;
 for(mt=1;mt<=ep;mt++)
      {s0=n/a;
        s1=0;
        for(k=0;k<b;k++)
             {i=k;
              while(i<n)
                	  {arg=i+b;
                		if(k==0)
                    	  {pr=xr[arg];
                    		pi=xi[arg];
            			  }
            			 else{
                                pr=xr[arg]*c[s1]-xi[arg]*s[s1];
                                pi=xr[arg]*s[s1]+xi[arg]*c[s1];
            					}
                         xr[arg]=xr[i]-pr;
                         xi[arg]=xi[i]-pi;
                         xr[i]=xr[i]+pr;
                         xi[i]=xi[i]+pi;
                		 i=i+a;
        			  }
              s1=s1+s0;
        	 }
      a=2*a;
      b=b*2;
      }
    free(c);
    free(s);
}
void tra(float *x,float *y)
{ float t;
  t=(*x);
  (*x)=(*y);
  (*y)=t;
}
void signal1(int n,float *yr,float *yi)
            	 {int i;
                    for(i=0;i<n;i++)
                         yr[i]=yi[i]=0;
                    for(i=0;i<4;i++)
                		 yr[i]=1;
        		 }
void signal2(int n,float *yr,float *yi)
            	{int i;
                for(i=0;i<n;i++)
                     yr[i]=yi[i]=0;
                for(i=0;i<4;i++)
                	 yr[i]=i+1;
                for(i=4;i<8;i++)
                	 yr[i]=8-i;
        		}
void signal3(int n,float *yr,float *yi)
            	{int i;
                 for(i=0;i<n;i++)
                      yr[i]=yi[i]=0;
                 for(i=0;i<4;i++)
                      yr[i]=4-i;
                 for(i=4;i<8;i++)
                      yr[i]=i-3;
        		}
void signal4(int n,float *yr,float *yi)
      {int i;
        for(i=0;i<n;i++)
             yr[i]=yi[i]=0;
        for(i=0;i<n;i++)
             yr[i]=cos(M_PI*i/4);
      }
void signal5(int n,float *yr,float *yi)
      {int i;
        for(i=0;i<n;i++)
            yr[i]=yi[i]=0;
        for(i=0;i<n;i++)
             yr[i]=sin(M_PI*i/8);
      }
void signal6(int n,float *yr,float *yi)
      {int i;
        for(i=0;i<n;i++)
             yr[i]=yi[i]=0;
        for(i=0;i<n;i++)
             yr[i]=0.2*cos(M_PI*i/8)+cos(M_PI*i/4)+cos(M_PI*i*20/64);
      }
void signal7(int n,float *yr,float *yi)
      {int i;
        for(i=0;i<n;i++)
             yr[i]=yi[i]=0;
        for(i=0;i<n;i++)
             yr[i]=cos(M_PI/4*i)+sin(M_PI/8*i);
      }
void signal8(int n,float *yr,float *yi)
      {int i;
        for(i=0;i<n;i++)
             yr[i]=cos(M_PI/4*i);
        for(i=0;i<n;i++)
             yi[i]=sin(M_PI/8*i);
      }
void main(void)
{ float y[64],yr[64],yi[64];
  int i,k=64;
  int m,n;
  InitGraphic();
  printf("ENTER SIGNAL(1-8):");
  scanf("%d",&m);
  if(m!=6)
     {printf("ENTER N:(8or16)");
      scanf("%d",&n);
     }
  else
  {printf("ENTER N:(16or32or64)");
    scanf("%d",&n);
  }
switch(m)
{case 1:signal1(n,yr,yi);break;
case 2:signal2(n,yr,yi);break;
case 3:signal3(n,yr,yi);break;
case 4:signal4(n,yr,yi);break;
case 5:signal5(n,yr,yi);break;
case 6:signal6(n,yr,yi);break;
case 7:signal7(n,yr,yi);break;
case 8:signal8(n,yr,yi);break;  }
PlotTimeSequ(yr,n,60,10,130);
FFT(yr,yi,n,0);
for(k=0;k<=n;k++) y[k]=sqrt(yr[k]*yr[k]+yi[k]*yi[k]);
PlotSpectrum(y,n,100,10,400);

setcolor(6);
outtextxy(200,440,"Press any key---");
getch();
closegraph();
return;
}