fft.c

一个基于QT的电网监控程序

/*
 * 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);
}