5.cpp

用最简单的程序去提取信号的功率谱

#include<stdlib.h>
#include<iostream.h>
#include<math.h>
#include<fstream.h>

void fft(double xreal[],double ximag[],int nfft,int sign);
void main()
{
	int i,n;
	double x[128],y[128];
	n=128;
	for(i=0;i<n;i++)
	{
      x[i]=1.0;
	  y[i]=0.0;
	}
	
	cout<<"-----1-----"<<endl;
	fft(x,y,n,1);
	cout<<"-----2-----"<<endl;
    ofstream file1;
	file1.open("a.txt",ios::out);
	for(i=0;i<128;i++)
	{  file1<<i<<"   "<<x[i]<<endl;
	    //cout<<i<<'\t'<<x[i]<<endl;
	}
	file1.close();
	cout<<"sucess!"<<endl;
}

void fft(double xreal[],double ximag[],int nfft,int sign)
{   
   int i,j,k,l,m,n,p,q;
   double tmp,tmp1,tmp2,tmp3,tmp4,nfft_2,e,pi,tr,ti;
   double c,s,w1,w2;
   pi=3.14159265;
   cout<<"------3-----"<<endl;
   nfft_2=log(nfft)/log(2);
   n=int(nfft_2);
   if(n<nfft_2)
     {
       n+=1;
       l=pow(2,n);
       for(i=nfft;i<l;i++)
          {xreal[i]=0.0;ximag[i]=0.0;}
       nfft=l;
      } ;
  /* for(i=0;i<nfft;i++)
	{      
	  cout<<"i="<<i<<'\t'<<xreal[i]<<'\t'<<ximag[i]<<endl;
	}*/
   for(i=0;i<nfft/4;i++)
     {	  
        k=2*i+1;
        j=k+(nfft/2-1);
		cout<<"k="<<k<<'\t'<<"j="<<j<<endl;
        tmp1=xreal[k];
        tmp2=ximag[k];
        xreal[k]=xreal[j];
        ximag[k]=ximag[j];
        xreal[j]=tmp1;
     }
   for(i=0;i<nfft;i++)
        ximag[j]=tmp2;
	{      
	  //cout<<"i="<<i<<'\t'<<xreal[i]<<'\t'<<ximag[i]<<endl;
	}
   //cout<<"------4-----"<<endl;
   for(l=0;l<n;l++)
      {  m=pow(2,l);
         cout<<"l="<<l<<endl;
		 cout<<"m="<<m<<endl;
		 p=m*2;
            for(i=0;i<nfft;i+=p)
                 {
				    q=i+m;
					cout<<"i="<<i<<endl;
					cout<<"m="<<m<<endl;
                    for(j=i;j<q;j++)
                       {
                          e=-sign*pi/m; 
                          c=cos(e);
                          s=sin(e);
                          w1=1.0;
                          w2=0.0;  
                          k=0;
                          while(k<j)
                               {
                                 tmp=w1;
                                 w1=w1*c-w2*s;
                                 w2=tmp*s+w2*c;
                                 k++;
								 //cout<<"w1="<<w1<<'\t'<<"w2="<<w2<<endl;
                                }
                          tmp1=xreal[j];
                          tmp2=ximag[j];
                          tmp3=xreal[j+m];
                          tmp4=ximag[j+m];
						  tr=w1*tmp3-w2*tmp4;
						  ti=w2*tmp3+w1*tmp4;
						  //cout<<tmp1<<'\t'<<tmp2<<'\t'<<tmp3<<'\t'<<tmp4<<endl;
                          xreal[j]=tmp1+(w1*tmp3-w2*tmp4);
                          ximag[j]=tmp2+(w2*tmp3+w1*tmp4);
                          xreal[j+m]=tmp1-(w1*tmp3-w2*tmp4);
                          ximag[j+m]=tmp2-(w2*tmp3+w1*tmp4);
						  cout<<j<<endl;
                        }
                  }
			cout<<"l++="<<l+1<<endl;
			cout<<"m="<<m<<endl;
        }
   cout<<"fft routing sucess!"<<endl;
     if(sign=-1)
       {
          for(i=0;i<nfft;i++)
             {
                xreal[i]=1.0/nfft*xreal[i];
                ximag[i]=1.0/nfft*ximag[i];
              }
        }
	 for(i=0;i<nfft;i++)
	 {   //cout<<"i="<<i<<'\t'<<xreal[i]<<'\t'<<ximag[i]<<endl;
	 }
 }