dsp.hpp

此程序基于c语言实现fft快速傅里叶变换

#include "Sequence.hpp"
#include "functions.hpp"
using namespace seq;
//在类定义外，要指明成员与友元函数
//--------------零化函数
//绝对值小于eps的都化为零 ,以后升级为美化函数! 
Sequence Sequence::zeros()
{
   int l=-left.size();
   int r=right.size();
   Sequence se=*this;
   for(int i=l;i<r;i++)
      if( abs(real( se[i] )) <eps )
         real(se[i])=0;
      else if( abs(imag(se[i])<eps) )
         imag(se[i])=0;
   return se;
}
         
//----------------卷积
//线性卷积  
Sequence seq::lincon(Sequence x,Sequence y)
{
    if( x.cycle()<=y.cycle() )
    {
        int lsize=x.leftsize()+y.leftsize();
        int rsize=x.rightsize()+y.rightsize()-1;
        Sequence linse(lsize,rsize);
        for(int n=-lsize;n<rsize;n++)
        {
            TP sum=0;
            for(int xi=-x.leftsize();xi<x.rightsize();xi++)
            {
                int yi=n-xi;
                if( yi>=y.leftsize() && yi<y.rightsize() )
                sum+=x[xi]*y[yi];
            }
            linse[n]=sum;
        }
        return linse.clean();
    }
    else return lincon(y,x);
}
//循环卷积 
Sequence seq::cyccon(Sequence x,Sequence y,int N)
{
   if( x.leftsize()+y.leftsize() )
   {
       cout <<"只有右边序列才能计算循环卷积!\n";
       system("pause");
       exit(1);
   }
   if( max(x.rightsize(),y.rightsize())>N )
   {
       cout <<"循环卷积长度应大于两序列的长度!\n";
       system("pause");
       exit(1);
   }
   Sequence cycse(0,N);
   for(int n=0;n<N;n++)
   {
       TP sum=0;
       for(int xi=0; xi<min(N,x.rightsize()); xi++)
       {
           int yi=(N+n-xi)%N;
           if( yi<min(N,y.rightsize()) )
           sum+=x[xi]*y[yi];
       }
       cycse[n]=sum;
   }
   return cycse.clean();
}
//-------------------------离散傅立叶变换
//------dft,按定义
Sequence seq::dft(Sequence x,int N=-1)
{
    if(x.leftsize()){
       cout <<"只计算因果序列的DFT!\n";
       system("pause");
       exit(1);
       }
    int n=x.rightsize();
    if(N<n && N>0){
       cout <<"DFT长度不能小于原序列长度!\n";
       system("pause");
       exit(1);
       }
    if(N>n) x[N-1];
    if(N==-1) N=x.rightsize();
    Sequence X(0,x.rightsize());
    for(int i=0;i<N;i++)
    {
       TP sum=0;
       for(int xi=0;xi<N;xi++)
          sum+=x[xi]*exp( TP(0,1)*TP(-1,0)*TP(2,0)*TP(3.14159,0)
                         /TP(N,0)*TP(i,0)*TP(xi,0) );
       X[i]=sum;
    }
    return X;
}
//------------fft & ifft , 2-dit
Sequence seq::ft(Sequence x,int N,string str)
{
   //输入参数检查 
   if(x.leftsize()){
       cout <<"只计算因果序列的FT!\n";
       system("pause");
       exit(1);
       }
   int n=x.rightsize();
   if(N<n && N>0){
       cout <<"FT长度不能小于原序列长度!\n";
       system("pause");
       exit(1);
       }
    //输入参数预处理   
   int power=0;
   if(N==0) N=x.rightsize();
   n=N; 
   n--;
   while(n)
   {  
      n>>=1;
      power++;
   }
   n=(1<<power);             //n为FFT序列长度 

   Sequence X[2];bool bX=0;
   X[0][n-1];X[1][n-1];
   for(int i=0;i<n;i++)
      X[0][i]=x[bitinv(power,i)];
   //FFT开始啦!
   //          蝶群数、    蝶群间距、每个蝶群碟数=每碟翅间距 
   unsigned butterflies, butterflies_d,    butterfly;
   butterflies=n>>1;
   butterflies_d=2;
   butterfly=1;
   while(butterflies)
   {  //蝶群循环 
      for(int bf=0;bf<butterflies;bf++)
      {  //蝶循环 
         for(int bfd=0;bfd<butterfly;bfd++)
         {  
            int xi=bf*butterflies_d+bfd;
            TP u,v;
            if(str=="fft"){
                 u=X[bX][xi];
                 v=X[bX][xi+butterfly]
                   *exp( TP(-1,0)*TP(0,1)*TP(2,0)*TP(3.14159,0)
                         /TP(butterflies_d,0)*TP(bfd,0) );
                 }
            else if(str=="ifft"){
                 u=X[bX][xi]/TP(2,0);
                 v=X[bX][xi+butterfly]
                   *exp( TP(1,0)*TP(0,1)*TP(2,0)*TP(3.14159,0)
                         /TP(butterflies_d,0)*TP(bfd,0)
                        )/TP(2,0);
                 }
            else {
                 cout <<"ft()调用错误\n";
                 system("pause");
                 exit(1);
                 }
            X[!bX][xi]=u+v;
            X[!bX][xi+butterfly]=u-v;
         }
      }
      butterflies>>=1;
      butterflies_d<<=1;
      butterfly<<=1;
      bX=!bX;
   }
   return X[bX];
}
//------------fft,2-dit-fft
Sequence seq::fft(Sequence x,int N=0)
{
 return ft(x,N,"fft");
}
//------------ifft,2-dit-ifft
Sequence seq::ifft(Sequence x,int N=0)
{
 return ft(x,N,"ifft");
}