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
