fft.cpp

fft代码,采用蝶形算法,包括C,matlab和verilog代码

// fft -- a program to model a Fast Fourier Transform (FFT).
//
// Copyright (C) 2001 John Dalton
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.



// Include required libraries

#include <vector>
#include <math.h>
#include <iostream>
#include <complex>
#include <stdlib.h>


//
// factor
//
// Return a factor of a number.
//
// Inputs:  x = the number to be factored.
// Returns: A factor of 'x'.  If 'x' is a
//          prime, eith 'x' or the number
//          one will be returned.
//
long factor(long x)
{
  long i;

  for(i=2; i<x; i++) {
    if(x%i==0) {
      return(i);
    }
  }
  return(1);
}


//Paramatise the type of each element
//to allow it to be changed easily

#define element_type double_complex


//This routine does one stage of an FFT.

//
// fft_stage
//
// Do a set of FFTs.  The length of the FFT may be
// smaller than the length of the input sequence.
// In this case, multiple FFTs will be performed.
// If the length of the FFT is a composite number,
// it is recusively decomposed into smaller FFTs
// whose lengths are prime numbers.  The elements
// of each FFT do not need to be consecutive.
//
// Inputs:  in = A sequence of numbers to be transformed.
//          fft_length = The number of points in each FFT
//          grouping = The separation between each element
//                     within an FFT.
// Returns: The transformed sequence.
//
vector<element_type> fft_stage(vector<element_type> in, int fft_length, int grouping)
{
  int x, y;
  vector<element_type> out(in.size());
  int f;

  //  cout<<"***FFT"<<fft_length<<" Grouping:"<<grouping<<" Array Length:"<<in.size()<<"\n";

  // A DFT of length does nothing, so just
  // return the input
  if(fft_length <=1) {
    out = in;
    return(out);
  }

  // Factorise the number of points in the DFT (if possible)
  // so we can break it into two smaller DFTs
  f = factor(fft_length);

  // If the number of points in the DFT is a composite
  // number (ie. can be factorised), divide it into
  // two smaller DFTs.  Recurse until the DFT has been
  // decomposed into DTFs of prime length.
  if(f!=fft_length &&f!=1) {
    int P = f;
    int Q = fft_length / f;
    vector<element_type> intermediate(in.size());

    //Do a DFT along one dimension
    intermediate = fft_stage(in, P, grouping*Q);
    //Multiply by twiddle factors
    for(int xx=0; xx<intermediate.size(); xx+=P*Q) {
      for(int a=0; a<P; a++) {
        for(int q=0; q<Q; q++) {
	  intermediate[xx+a*Q+q] = intermediate[xx+a*Q+q] * exp(-double_complex(0,1)*2*M_PI/(P*Q)*a*q);
        }
      }
    }
    //Do a DFT along a second, orthogonal, dimension
    out = fft_stage(intermediate, Q, grouping);

    return(out);
  }

  // Here if the length for the DFT was prime.   
  for(x=0; x<in.size(); x+=fft_length*grouping) {
    for(y=x; y<x+grouping; y++) {
      //If the length of the DFT is two, do a butterfly
      //operation.
      if(fft_length==2) {  
        out[y         ] = in[y] + in[y+grouping]; 
        out[y+grouping] = in[y] - in[y+grouping];
      } else { //do an exhaustive DFT (for the time being)
        //For other prime length DTFs (3,5,7,11,...) do a
        //full DFT.  I'm still looking for a more efficient
        //algorithm to use here.  Any suggestions?
        for(int k=0; k<fft_length; k++) {
          out[y+k*grouping] = 0;
          for(int n=0; n<fft_length; n++){
            out[y+k*grouping] += in[y+n*grouping]*exp(-double_complex(0,1)*n*k*2*M_PI/fft_length);
	  }
	}
      }
    }
  }
  return(out);
}



//
// fft
//
// Do an FFT.  There are no restrictions on the
// length of the FFT, but a power of two will
// yield the most efficient calculation.
// This is the routine you should call to do
// an FFT (not fft_stage).
//
// Inputs:  in = A sequence of numbers to be transformed.
// Returns: The transformed sequence.
//
vector<element_type> fft(vector<element_type> in)
{
  return( fft_stage(in, in.size(), 1) );
}



//
// golden_dft_model
//
// Do a Discrete Fourier Transform (DFT).
// This is the model against which results
// of the FFT will be compared.
//
// Inputs:  in = A sequence of numbers to be transformed.
// Returns: The transformed sequence.
//
vector<element_type> golden_dft_model(vector<element_type> in)
{
  vector<element_type> out(in.size());

  for(int k=0; k<out.size(); k++) {
    out[k] = 0;
    for(int n=0; n<in.size(); n++) {
      out[k] += in[n]*exp(-double_complex(0,1)*2*M_PI*n*k/in.size());
    }
  }

  return(out);
}


//Reproduce the order in which the outputs of the FFT appear
//
// order
//
// This function reproduces the order of the
// outputs of the FFT.  For an FFT of length a
// power of two, the order will be bit reversed.
// For other lengths, the ordering will be slightly
// more complicated.
//
// Inputs:  in = A sequence of numbers 1,2,3,...,N-1.
//          length = The length of the FFT whose output
//                   ordering is to be reproduced.
//          grouping = The separation between each element
//                     of the FFT.
// Returns: A sequence whose order corresponds to the
//          order of the FFT's outputs.
//
vector<int> order(vector<int> in, int length, int grouping)
{
  vector<int> out(in.size());
  int f;
 
  if(length <=1) {
    out = in;
    return(out);
  }

  f = factor(length);

  if(f!=length &&f!=1) {
    int P = f;
    int Q = length / f;
    vector<int> intermediate(in.size());
    //    intermediate = order(in, P, grouping*Q);
    //out = order(intermediate, Q, grouping);
    intermediate = order(in, Q, grouping*P);
    out = order(intermediate, P, grouping);
    return(out);
  }

  for(int x=0; x<in.size(); x+=length*grouping) {
    for(int y=0; y<grouping; y++) {
      for(int k=0; k<length; k++) {
        out[x+y+k*grouping] = in[x+y*length+k];
      }
    }
  }
  return(out);
}



//
// near
//
// Check if two complex numbers are equal,
// to within some approximation.
//
// Inputs:  a = The first complex number.
//          b = The second complex number.
// Returns: One if the two numbers are equal, zero otherwise.
//
int near(element_type a, element_type b)
{
  const double epsilon = 1e-12;
  if( fabs(a.real()-b.real()) <= epsilon && fabs(a.imag()-b.imag()) <= epsilon ) {
    return(1);
  }
  return(0);
}


//
// test_fft
//
// Do an FFT and compare the results with the
// output of the corresponging DFT.  If the
// results are not the same, print an
// error message.
//
// Inputs:  in = The complex sequence to be transformed.
// Outputs: failed = Non-zero if the test failed, zero if it suceeded.
// Returns: The results of the FFT.
//
vector<element_type> test_fft(vector<element_type> in, int *failed)
{
  vector<element_type> out(in.size()), reference(in.size());
  vector<int> index(in.size());

  out = fft(in);
  reference = golden_dft_model(in);

  for(int i=0; i<index.size(); i++) {
    index[i] = i;
  }

  index = order(index, index.size(), 1);
   
  //for(int i=0; i<index.size(); i++) {
  //  cout << index[i] << ((i==index.size()-1)?'\n':',');
  //}

  for(int i=0; i<index.size(); i++) {
    if( !near(reference[i], out[index[i]]) ) {
      cout << "FFT Test Failed!\n";
      cout << i << ": " << reference[i] << " != " << out[index[i]] << "\n";
      for(int j=0; j<in.size(); j++) {
        cout << j << ": " << in[j] << " " << reference[j] << " " << out[index[j]] << "\n";
      }
      *failed = 1;
      return(out);
    }
  }

  *failed = 0;
  return(out);
}


//
// main
//
// Run a series of tests to verify an FFT.
//
// Inputs:  None.
// Outputs: None.
// Returns: Zero if the tests suceed.  Non-zero otherwise.
//
main()
{
  int i;
  int length;
  int failed;
  //const double freq = 3.0;

  // Do a series of tests using walking tones on a
  // series of different length FFTs.  Assuming linearity,
  // these tests should prove the FFT correct.
  cout<< "Tone Test.\n";
  cout.flush();
  for(length=1; length<=32; length++) {
    vector<element_type> in(length), out(length);
    for(double freq=0; freq<length; freq++){

      for(i=0; i<length; i++) {
        //  in[i] = i;
        in[i] = exp(double_complex(0,1)*2.0*M_PI/length*freq*i);
      }

      for(i=0; i<length; i++) {
        //   cout << in[i] << ((i==length-1)?'\n':',');
      }

      out = test_fft(in, &failed);
      if(failed) {
        exit(1);
      }

      for(i=0; i<length; i++) {
        //cout << "(" << (int)(out[i].real()*100.0)/100 << "," << (int)(out[i].imag()*100.0)/100 << ")"  << ((i==length-1)?'\n':',');
      }
    }
    }

  // Do a series of tests with random data and
  // random length FFTs.  These test double check
  // that nothing is wrong.
  cout<< "Monte Carlo Test.\n";
  cout.flush();
  int max_length = 100;
  for(i=0; i<100000; i++) {
    int length = rand()%max_length;
    vector<element_type> in(length);
    //   cout<< length << " ";
    //cout.flush();
    for(int i=0; i<in.size(); i++) {
      in[i] = 1.0-2.0*(double)rand()/(double)RAND_MAX;
    }
    test_fft(in, &failed);
    if(failed) {
      exit(1);
    }
  }

  cout << "PASSED ALL TESTS.\n";
  return(0);
}