//
// Programmer:    Craig Stuart Sapp <craig@ccrma.stanford.edu>
// Creation Date: Fri Feb  7 00:44:33 PST 1997
// Last Modified: Thu Dec 18 18:39:40 GMT-0800 1997
// Last Modified: Sun Jun 11 12:55:07 PDT 2000 (improved memory management)
// Last Modified: Sun Jul  8 14:58:52 PDT 2001 (removed templates)
// Filename:      ...sig/maint/code/utilities/transforms/transforms.cpp
// Web Address:   http://sig.sapp.org/src/sig/transforms.cpp
// Syntax:        C++ 
//
// Description:   Basic signal <-> frequency transforms.
//

#ifndef _TRANSFORMS_CPP_INCLUDED
#define _TRANSFORMS_CPP_INCLUDED

#include "Block.h"
#include "ComplexD.h"

#include "transforms.h"
#include "transforms-private.h"

#include <math.h>
#include <stdlib.h>
#include <string.h>

#ifndef OLDCPP
   #include <iostream>
   using namespace std;
#else
   #include <iostream.h>
#endif


//////////////////////////////
//
// FFT -- Fast Fourier Transform O(N Log N)
//   Returns the complex spectrum of the given complex input signal.
//   Length of Block must be a power of 2.
//

Block<ComplexD>& FFT(Block<ComplexD>& outputSpectrum, 
      const Block<ComplexD>& inputSignal) {
   int N = inputSignal.getSize();
 
   if (!poweroftwo(N)) {
      cerr << "You can only take the FFT of a block with length being" 
           << " a power of 2.\nRequested transform length: " << N << endl;
      exit(1);
   }

   outputSpectrum = inputSignal;
   fft_destructive(outputSpectrum);
   return outputSpectrum;
}



//////////////////////////////
//
// IFFT -- Inverse Fast Fourier Transform O(N Log N).
//   Returns the complex signal of the given complex input spectrum.
//   Length of Block must be a power of 2.
//

Block<ComplexD>& IFFT(Block<ComplexD>& outputSignal, 
      const Block<ComplexD>& inputSpectrum) {
   int N = inputSpectrum.getSize();
 
   if (!poweroftwo(N)) {
      cerr << "You can only take the IFFT of a block with length being" 
           << " a power of 2.\nRequested transform length: " << N << endl;
      exit(1);
   }

   outputSignal = inputSpectrum;
   ifft_destructive(outputSignal);
   return outputSignal;
}



//////////////////////////////
//
// DirectDFT -- Direct Discrete Fourier Transform O(N^2).
//   Returns the complex spectrum of the given complex input signal.
//   Quickly becomes slower than the FFT if the length of the input
//   block is greater than about 100.
//

Block<ComplexD>& DirectDFT(Block<ComplexD>& outputSpectrum, 
      const Block<ComplexD>& inputSignal) {
   int N = inputSignal.getSize();
   const Block<ComplexD>& x = inputSignal;
   outputSpectrum.setSize(N);
   Block<ComplexD>& X = outputSpectrum;
   ComplexD sum;
   int k, n;
   ComplexD kernel(0,0);
   double mtpn = -2 * 3.14159265359 / N;


   for (k=0; k<N; k++) {
      sum.zero();
      for (n=0; n<N; n++) {
        kernel.im() = mtpn * k * n;
        sum += x[n] * kernel.cexp();
      }
      X[k] = sum;
   }
  
   return outputSpectrum;
}



////////////////////////////// 
//
// DFT -- Selects the DirectDFT or the FFT as appropriate.
//

Block<ComplexD>& DFT(Block<ComplexD>& outputSpectrum, 
      const Block<ComplexD>& inputSignal) {
   int N = inputSignal.getSize();
   if (poweroftwo(N) && N > 63) {
      return FFT(outputSpectrum, inputSignal);
   } else {
      return DirectDFT(outputSpectrum, inputSignal);
   }
}



//////////////////////////////
//
// DirectIDFT -- Direct Inverse Discrete Fourier Transform O(N^2).
//   Returns the complex spectrum of the given complex input signal.
//   Quickly becomes slower than the IFFT if the length of the input
//   block is greater than about 100.
//

Block<ComplexD>& DirectIDFT(Block<ComplexD>& outputSignal, 
      const Block<ComplexD>& inputSpectrum) {
   int N = inputSpectrum.getSize();
   const Block<ComplexD>& X = inputSpectrum;
   outputSignal.setSize(N);
   Block<ComplexD>& x = outputSignal;
   ComplexD sum;
   int k, n;
   ComplexD kernel(0,0);
   double tpn = 2.0 * 3.14159265359 / (double)N;

   for (n=0; n<N; n++) {
      sum.zero();
      for (k=0; k<N; k++) {
        kernel.im() = tpn * k * n;
        sum += X[k] * kernel.cexp();
      }
      x[n] = sum/N;
   }
  
   return outputSignal;
}



////////////////////////////// 
//
// IDFT -- Selects the DirectIDFT or the IFFT as appropriate.
//

Block<ComplexD>& IDFT(Block<ComplexD>& outputSignal, 
      const Block<ComplexD>& inputSpectrum) {
   int N = inputSpectrum.getSize();
   if (poweroftwo(N) && N > 63) {
      return IFFT(outputSignal, inputSpectrum);
   } else {
      return DirectIDFT(outputSignal, inputSpectrum);
   }
}



#endif  /* _TRANSFORMS_CPP_INCLUDED */



// md5sum: 14d6cb9182c5c4eca979b24a829214b7 transforms.cpp [20050403]