/* * Xavier Gourdon : Sept. 99 (xavier.gourdon.free.fr) *  * FFT.c : Very basic FFT file. *         The FFT encoded in this file is short and very basic. *         It is just here as a teaching file to have a first  *         understanding of the FFT technique. *          * A lot of optimizations could be made to save a factor 2 or 4 for time * and space. *  - Use a 4-Step (or more) FFT to avoid data cache misses. *  - Use an hermitian FFT to take into account the hermitian property of *    the FFT of a real array. *  - Use a quad FFT (recursion N/4->N instead of N/2->N) to save 10 or *    15% of the time. * *  Informations can be found on  *    http://xavier.gourdon.free.fr/Constants/constants.html */ #include <stdlib.h>#include <stdio.h>#include <math.h>#include "fft.h"long FFTLengthMax;Complex * OmegaFFT;Complex * ComplexCoef; double FFTSquareWorstError;long AllocatedMemory;/*初始化*/void InitializeFFT(long MaxLength){  long i;  Real Step;  FFTLengthMax = MaxLength;  OmegaFFT = (Complex *) malloc(MaxLength/2*sizeof(Complex));  ComplexCoef = (Complex *) malloc(MaxLength*sizeof(Complex));  Step = 2.*PI/(double) MaxLength;  for (i=0; 2*i<MaxLength; i++) {    OmegaFFT[i].R = cos(Step*(double)i);    OmegaFFT[i].I = sin(Step*(double)i);  }  FFTSquareWorstError=0.;  AllocatedMemory = 7*MaxLength*sizeof(Complex)/2;}/*遞歸FFT*/void RecursiveFFT(Complex * Coef, Complex * FFT, long Length, long Step, long Sign){  long i, OmegaStep;  Complex * FFT0, * FFT1, * Omega;  Real tmpR, tmpI;  if (Length==2) {    FFT[0].R = Coef[0].R + Coef[Step].R;    FFT[0].I = Coef[0].I + Coef[Step].I;    FFT[1].R = Coef[0].R - Coef[Step].R;    FFT[1].I = Coef[0].I - Coef[Step].I;    return;  }  FFT0 = FFT;  RecursiveFFT(Coef     ,FFT0,Length/2,Step*2,Sign);  FFT1 = FFT+Length/2;  RecursiveFFT(Coef+Step,FFT1,Length/2,Step*2,Sign);  Omega = OmegaFFT;  OmegaStep = FFTLengthMax/Length;  for (i=0; 2*i<Length; i++, Omega += OmegaStep) {    /* 遞歸公式:       FFT[i]          <-  FFT0[i] + Omega*FFT1[i]       FFT[i+Length/2] <-  FFT0[i] - Omega*FFT1[i],       Omega = exp(2*I*PI*i/Length) */    tmpR = Omega[0].R*FFT1[i].R-Sign*Omega[0].I*FFT1[i].I;    tmpI = Omega[0].R*FFT1[i].I+Sign*Omega[0].I*FFT1[i].R;    FFT1[i].R = FFT0[i].R - tmpR;    FFT1[i].I = FFT0[i].I - tmpI;    FFT0[i].R = FFT0[i].R + tmpR;    FFT0[i].I = FFT0[i].I + tmpI;  }}/* 計算 Coef 的復數傅立葉變換 FFT */void FFT(Real * Coef, long Length, Complex * FFT, long NFFT){  long i;  /* 實數系數傳給復數系數的實部 */  for (i=0; i<Length; i++) {    ComplexCoef[i].R = Coef[i];    ComplexCoef[i].I = 0.;  }  for (; i<NFFT; i++)    ComplexCoef[i].R = ComplexCoef[i].I = 0.;  RecursiveFFT(ComplexCoef,FFT,NFFT,1,1);}void EndFFT(){ free(OmegaFFT); free(ComplexCoef);}