FFT算法的c语言实现

// 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 
#include 
#include 
#include 
#include 


//
// 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    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 fft_stage（vector in int fft_length int grouping）
{
  int x y;
  vector out（in.size（））;
  int f;

  //  cout<<“***FFT“<
  // 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 intermediate（in.size（））;

    //Do a DFT along one dimension
    intermediate = fft_stage（in P grouping*Q）;
    //Multiply by twiddle fact