libs/qm-dsp/dsp/transforms/FFT.cpp

   1 /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
   2
   3 /*
   4     QM DSP Library
   5
   6     Centre for Digital Music, Queen Mary, University of London.
   7     This file is based on Don Cross's public domain FFT implementation.
   8 */
   9
  10 #include "FFT.h"
  11
  12 #include "maths/MathUtilities.h"
  13
  14 #include <cmath>
  15
  16 #include <iostream>
  17
  18 FFT::FFT(unsigned int n) :
  19     m_n(n),
  20     m_private(0)
  21 {
  22     if( !MathUtilities::isPowerOfTwo(m_n) )
  23     {
  24         std::cerr << "ERROR: FFT: Non-power-of-two FFT size "
  25                   << m_n << " not supported in this implementation"
  26                   << std::endl;
  27         return;
  28     }
  29 }
  30
  31 FFT::~FFT()
  32 {
  33
  34 }
  35
  36 FFTReal::FFTReal(unsigned int n) :
  37     m_n(n),
  38     m_private_real(0)
  39 {
  40     m_private_real = new FFT(m_n);
  41 }
  42
  43 FFTReal::~FFTReal()
  44 {
  45     delete (FFT *)m_private_real;
  46 }
  47
  48 void
  49 FFTReal::process(bool inverse,
  50                  const double *realIn,
  51                  double *realOut, double *imagOut)
  52 {
  53     ((FFT *)m_private_real)->process(inverse, realIn, 0, realOut, imagOut);
  54 }
  55
  56 static unsigned int numberOfBitsNeeded(unsigned int p_nSamples)
  57 {
  58     int i;
  59
  60     if( p_nSamples < 2 )
  61     {
  62         return 0;
  63     }
  64
  65     for ( i=0; ; i++ )
  66     {
  67         if( p_nSamples & (1 << i) ) return i;
  68     }
  69 }
  70
  71 static unsigned int reverseBits(unsigned int p_nIndex, unsigned int p_nBits)
  72 {
  73     unsigned int i, rev;
  74
  75     for(i=rev=0; i < p_nBits; i++)
  76     {
  77         rev = (rev << 1) | (p_nIndex & 1);
  78         p_nIndex >>= 1;
  79     }
  80
  81     return rev;
  82 }
  83
  84 void
  85 FFT::process(bool p_bInverseTransform,
  86              const double *p_lpRealIn, const double *p_lpImagIn,
  87              double *p_lpRealOut, double *p_lpImagOut)
  88 {
  89     if (!p_lpRealIn || !p_lpRealOut || !p_lpImagOut) return;
  90
  91 //    std::cerr << "FFT::process(" << m_n << "," << p_bInverseTransform << ")" << std::endl;
  92
  93     unsigned int NumBits;
  94     unsigned int i, j, k, n;
  95     unsigned int BlockSize, BlockEnd;
  96
  97     double angle_numerator = 2.0 * M_PI;
  98     double tr, ti;
  99
 100     if( !MathUtilities::isPowerOfTwo(m_n) )
 101     {
 102         std::cerr << "ERROR: FFT::process: Non-power-of-two FFT size "
 103                   << m_n << " not supported in this implementation"
 104                   << std::endl;
 105         return;
 106     }
 107
 108     if( p_bInverseTransform ) angle_numerator = -angle_numerator;
 109
 110     NumBits = numberOfBitsNeeded ( m_n );
 111
 112
 113     for( i=0; i < m_n; i++ )
 114     {
 115         j = reverseBits ( i, NumBits );
 116         p_lpRealOut[j] = p_lpRealIn[i];
 117         p_lpImagOut[j] = (p_lpImagIn == 0) ? 0.0 : p_lpImagIn[i];
 118     }
 119
 120
 121     BlockEnd = 1;
 122     for( BlockSize = 2; BlockSize <= m_n; BlockSize <<= 1 )
 123     {
 124         double delta_angle = angle_numerator / (double)BlockSize;
 125         double sm2 = -sin ( -2 * delta_angle );
 126         double sm1 = -sin ( -delta_angle );
 127         double cm2 = cos ( -2 * delta_angle );
 128         double cm1 = cos ( -delta_angle );
 129         double w = 2 * cm1;
 130         double ar[3], ai[3];
 131
 132         for( i=0; i < m_n; i += BlockSize )
 133         {
 134
 135             ar[2] = cm2;
 136             ar[1] = cm1;
 137
 138             ai[2] = sm2;
 139             ai[1] = sm1;
 140
 141             for ( j=i, n=0; n < BlockEnd; j++, n++ )
 142             {
 143
 144                 ar[0] = w*ar[1] - ar[2];
 145                 ar[2] = ar[1];
 146                 ar[1] = ar[0];
 147
 148                 ai[0] = w*ai[1] - ai[2];
 149                 ai[2] = ai[1];
 150                 ai[1] = ai[0];
 151
 152                 k = j + BlockEnd;
 153                 tr = ar[0]*p_lpRealOut[k] - ai[0]*p_lpImagOut[k];
 154                 ti = ar[0]*p_lpImagOut[k] + ai[0]*p_lpRealOut[k];
 155
 156                 p_lpRealOut[k] = p_lpRealOut[j] - tr;
 157                 p_lpImagOut[k] = p_lpImagOut[j] - ti;
 158
 159                 p_lpRealOut[j] += tr;
 160                 p_lpImagOut[j] += ti;
 161
 162             }
 163         }
 164
 165         BlockEnd = BlockSize;
 166
 167     }
 168
 169
 170     if( p_bInverseTransform )
 171     {
 172         double denom = (double)m_n;
 173
 174         for ( i=0; i < m_n; i++ )
 175         {
 176             p_lpRealOut[i] /= denom;
 177             p_lpImagOut[i] /= denom;
 178         }
 179     }
 180 }
 181