Implementation of several fast Fourier transform (FFT) C + + implementations

Source: Internet
Author: User
Tags cos sin

Links:http://blog.csdn.net/zwlforever/archive/2008/03/14/2183049.aspx
A good FFT article, a collection of a bit.
The positive and inverse transformations of DFT are (1) and (2) respectively. Assuming that there are n data, then the calculation of a frequency point requires N-times complex multiplication and N-1 complex addition, the entire DFT requires n*n complex multiplication and n (N-1) complex addition; Since the multiplication of complex numbers requires 4 real multiplication and 2 complex additions, the first complex addition takes two real additions, Therefore the whole DFT needs 4*n*n of the real number multiplication and 2*n (N-1) +2*n*n≈4*n*n the plural addition of the times. When n is large, the required computational effort is quite large, for example, it takes about 4 million multiplication when n=1024, for real-time signal processing, the computing equipment will be very demanding, so it is proposed how to reduce the computational DFT computation problem.

in 1965, the Library Force and EOG published in the journal Computational Mathematics an algorithm for calculating Fourier series of machines, which was first proposed by FFT algorithm. Since then, the fast algorithm of DFT is called hot topic of research, and it is also the advent of FFT algorithm, which makes digital signal processing can be applied to real-time situations and further promote the development and application of digital signal processing.

Most FFT algorithms use the Characteristics of (3) Periodic, conjugate symmetry, compressibility and extensibility to compress the computational amount.

1) code calculated according to the DfT definition

//data is the input data pointer, log2n=log2 (length), Flag=-1 is a positive transformation, flag=1 is the inverse transform, the transformation results are also returned by the pointer data point to the space
void DfT (complex<double>*data,int log2n,int flag)
{
int i,j,length;
Complex<double> WN;
length=1<<log2n;
Complex<double>*temp=new complex<double> (length);
for (i=0;i<length;i++)
{
temp[i]=0;
for (j=0;j<length;j++)
{
                     wn=complex<double> (cos (2.0*pi/length*i*j), sin (flag*2.0* PI/LENGTH*I*J));
                     temp[i]+=data[j]*wn;    
                 }            
       }
       if ( flag==1)
           for (i=0;i<length;i++)
               data[i]=temp[i]/ Length
        delete[] temp;
}   2), base 2 time selected FFT

The digital signal sequence of time domain is grouped by odd and even, the following transformations can be made, and from the transformation results we can know that a DFT of length N can be converted into a combination of two sub-sequences of length N/2. And so on, it is possible to convert a combination of Fourier changes to N/2 2 points. However, the input should be a 2-based inverted order.

As a result of a number of parity, input data to be in accordance with the principle of 2 reverse order, the output data is the normal sequence, the reverse algorithm is described in addition. The following first describes the algorithm in the form of recursion, because the recursive method does not fully utilize the advantages of the DIT2 algorithm---in-situ calculation, so the recursive form is only for the sake of clarity.

void Dit2rec (Complex<double>*indata,complex<double>*outdata,int length,int sign)
{
Complex<double>*evendata=new complex<double> (LENGTH/2);
Complex<double>*odddata =new complex<double> (LENGTH/2);
Complex<double>*evenresult=new complex<double> (LENGTH/2);
Complex<double>*oddresult=new complex<double> (LENGTH/2);
int i,j;
if (length==1)
{
      outdata[0]=indata[0]/n;
      return;
   }
  for (i=0;i<length/2;i++)
  {
     EvenData[i]=InData[2*i];
    OddData[i]=InData[2*i+1];
  }
  dit2rec (evendata,evenresult,length/2,sign);
  dit2rec (odddata,evenresult,length/2,sign);
  for (i=0;i<length/2;i++)
  {
     OutData[i]=EvenData+OddData*complex<double> (cos (2*pi*i/length), sin (sign*2*pi*i/length));
    OutData[i+length/2]=EvenData- OddData*complex<double> (cos (2*pi*i/length ), sin (sign*2*pi*i/length));
  }
  delete[] EvenData,OddData,EvenResult,OddResult;
  return;
}
 the non-recursive implementation is as follows (the reverse number of input is not considered):void Dit2 (Complex<double>*data,int log2n,int sign)
{
int i,j,k,step,length;
Complex<double> Wn,temp,deltawn;
length=1<<log2n;
for (i=1;i<=log2n;i++)
{
wn=1;step=1<<i;deltawn=complex<double> (cos (2*pi/step), sin (sign*2.0*pi/step));
for (j=0;j<step/2;j++)
{&NBSP;&NBSP;&NBSP;&NBSP;&NBSP;&NBSP;&NBSP;&NBSP;
         for (i=0;i<length/step;i++)
        {
           temp=data[i* Step+step/2+j]*wn;
           data[i*step+step/2+j]=data[i*step+j]- Temp
           Data[i*step+j]=Data[i*step+j]+temp;
         }
        wn=wn*deltawn;
      }
    }
    if (sign==-1)
       for (i=0;i<length;i++)
            Data[i]/=length;
 }

when I=1, that is, the first cycle is not necessary to do complex operations, because J can only take 1,wn as the real number, this time can be saved. This can therefore be improved to:void Dit2 (Complex<double>*data,int log2n,int sign)
{
int i,j,k,step,length;
Complex<double> Wn,temp,deltawn;
length=1<<log2n;
for (i=0;i<length;i+=2)
{
Temp=data[i];
DATA[I]=DATA[I]+DATA[I+1];
DATA[I+1]=TEMP-DATA[I+1];
}
for (i=2;i<=log2n;i++)
{
wn=1;step=1<<i;deltawn=complex<double> (cos (2.0*pi/step), sin (sign*2.0*pi/step));;
for (j=0;j<step/2;j++)
{&NBSP;&NBSP;&NBSP;&NBSP;&NBSP;&NBSP;&NBSP;&NBSP;
         for (i=0;i<length/step;i++)
        {
           temp=data[i* Step+step/2+j]*wn;
           data[i*step+step/2+j]=data[i*step+j]- Temp
           Data[i*step+j]=Data[i*step+j]+temp;
         }
         wn=wn*deltawn;
      }
   }
   if (sign==1)
   for (i=0;i<length;i++)
      Data[i]/=length;
}

3), base 2 frequency decimation FFTthe recursive version of the//DIF2 implementation:
void Dif2rec (Complex<double>*indata,complex<double>*outdata,int length,int sign)
{
complex<double>* ldata=new complex<double> (LENGTH/2);
complex<double>* lresult=new complex<double> (LENGTH/2);
complex<double>* rdata=new complex<double> (LENGTH/2);
complex<double>* rresult=new complex<double> (LENGTH/2);
Complex<double> temp;
int i;
if (length==1)
{
       outdata[0]=indata[0];
       return;
}
for (i=0;i<length/2;i++)
{
   ldata[i]=indata[i] ;
   RData[i]=InData[i+length/2];
}
for (i=0;i<length/2;i++)
{
   temp=ldata[i];
   LData[i]=LData[i]+RData[i];
   rdata[i]= (Temp-rdata[i]) * complex<double> (cos (2*pi*i/length), sin (sign*2*pi*i/ Length))
}
Dit2rec (ldata,lresult,length/2,sign);
Dit2rec (rdata,rresult,length/2,sign);
   for (i=0;i<length/2;i++)
   {
      OutData[2*i]=LResult[i];
      OutData[2*i+1]=RResult[i];
}
return;
}
  the//non-recursive implementation is as follows (the reverse number of input is not considered):
void Dif2 (Complex<double>*indata,int r,int sign)
{
Int length=1<<r;
int i,j,k,step;
complex<double> temp,wn;
for (i=1;i<=r;i++)
{
   step=1<< (r-i+1);
   wn=1;
   for (j=0;j<step/2;j++)
   {
       for (k=0;k<step/2;k++)
      {
         temp=indata[k*step+j];
         indata[k*step+j]=indata[k*step+j]-indata[k*step+step/2+j ];
         indata[k*step+step/2+j]= (temp-InData[k*step+step/2+j ]) *wn;
}
wn=wn*complex<double> (cos (2*pi/step*j), sin (sign*2*pi/step*j));
}
}
}
 

As with the DIT, the last loop of the outermost layer can be separate, because the last loop does not have to be complex, which reduces the number of complex operations.

A fast Fourier algorithm based on four-time decimation

Implementation of several fast Fourier transform (FFT) C + + implementations

Contact Us

The content source of this page is from Internet, which doesn't represent Alibaba Cloud's opinion; products and services mentioned on that page don't have any relationship with Alibaba Cloud. If the content of the page makes you feel confusing, please write us an email, we will handle the problem within 5 days after receiving your email.

If you find any instances of plagiarism from the community, please send an email to: info-contact@alibabacloud.com and provide relevant evidence. A staff member will contact you within 5 working days.

A Free Trial That Lets You Build Big!

Start building with 50+ products and up to 12 months usage for Elastic Compute Service

  • Sales Support

    1 on 1 presale consultation

  • After-Sales Support

    24/7 Technical Support 6 Free Tickets per Quarter Faster Response

  • Alibaba Cloud offers highly flexible support services tailored to meet your exact needs.