Code cleanup: removed the old real FFT code not used since at least 2009.

I confirmed that the currently used real FFT code in RealFFTf.cpp is faster
than the old one with a quick benchmark that calls PowerSpectrum() on 4-minute
audio file, with different sizes of computation windows:

Window_size: 256 method: new FFT time_s: 0.393
Window_size: 256 method: old FFT time_s: 1.065
Window_size: 1024 method: new FFT time_s: 0.38
Window_size: 1024 method: old FFT time_s: 0.958
Window_size: 4096 method: new FFT time_s: 0.413
Window_size: 4096 method: old FFT time_s: 1.084
Window_size: 16384 method: new FFT time_s: 0.518
Window_size: 16384 method: old FFT time_s: 1.338
Window_size: 65536 method: new FFT time_s: 0.655
Window_size: 65536 method: old FFT time_s: 1.524
Window_size: 262144 method: new FFT time_s: 0.735
Window_size: 262144 method: old FFT time_s: 1.873

Author: RaphaelMarinier

src/Experimental.h

@@ -107,10 +107,6 @@
 // Paul Licameli (PRL) 29 Nov 2014
 // #define EXPERIMENTAL_IMPROVED_SEEKING
 
-// Philip Van Baren 01 July 2009
-// Replace RealFFT() and PowerSpectrum function to use (faster) RealFFTf function
-#define EXPERIMENTAL_USE_REALFFTF
-
 // RBD, 1 Sep 2008
 // Enables MIDI Output of NoteTrack (MIDI) data during playback
 // USE_MIDI must be defined in order for EXPERIMENTAL_MIDI_OUT to work

src/FFT.cpp

@@ -46,6 +46,7 @@
 #include <stdio.h>
 #include <math.h>
 
+#include "RealFFTf.h"
 #include "Experimental.h"
 
 static int **gFFTBitTable = NULL;
@@ -110,10 +111,6 @@ void InitFFT()
    }
 }
 
-#ifdef EXPERIMENTAL_USE_REALFFTF
-#include "RealFFTf.h"
-#endif
-
 void DeinitFFT()
 {
    if (gFFTBitTable) {
@@ -122,10 +119,8 @@ void DeinitFFT()
       }
       delete[] gFFTBitTable;
    }
-#ifdef EXPERIMENTAL_USE_REALFFTF
    // Deallocate any unused RealFFTf tables
    CleanupFFT();
-#endif
 }
 
 inline int FastReverseBits(int i, int NumBits)
@@ -238,22 +233,11 @@ void FFT(int NumSamples,
 /*
  * Real Fast Fourier Transform
  *
- * This function was based on the code in Numerical Recipes in C.
- * In Num. Rec., the inner loop is based on a single 1-based array
- * of interleaved real and imaginary numbers.  Because we have two
- * separate zero-based arrays, our indices are quite different.
- * Here is the correspondence between Num. Rec. indices and our indices:
- *
- * i1  <->  real[i]
- * i2  <->  imag[i]
- * i3  <->  real[n/2-i]
- * i4  <->  imag[n/2-i]
+ * This is merely a wrapper of RealFFTf() from RealFFTf.h.
  */
 
 void RealFFT(int NumSamples, float *RealIn, float *RealOut, float *ImagOut)
 {
-#ifdef EXPERIMENTAL_USE_REALFFTF
-   // Remap to RealFFTf() function
    int i;
    HFFT hFFT = GetFFT(NumSamples);
    float *pFFT = new float[NumSamples];
@@ -280,62 +264,8 @@ void RealFFT(int NumSamples, float *RealIn, float *RealOut, float *ImagOut)
    }
    delete [] pFFT;
    ReleaseFFT(hFFT);
-
-#else
-
-   int Half = NumSamples / 2;
-   int i;
-
-   float theta = M_PI / Half;
-
-   float *tmpReal = new float[Half];
-   float *tmpImag = new float[Half];
-
-   for (i = 0; i < Half; i++) {
-      tmpReal[i] = RealIn[2 * i];
-      tmpImag[i] = RealIn[2 * i + 1];
-   }
-
-   FFT(Half, 0, tmpReal, tmpImag, RealOut, ImagOut);
-
-   float wtemp = float (sin(0.5 * theta));
-
-   float wpr = -2.0 * wtemp * wtemp;
-   float wpi = -1.0 * float (sin(theta));
-   float wr = 1.0 + wpr;
-   float wi = wpi;
-
-   int i3;
-
-   float h1r, h1i, h2r, h2i;
-
-   for (i = 1; i < Half / 2; i++) {
-
-      i3 = Half - i;
-
-      h1r = 0.5 * (RealOut[i] + RealOut[i3]);
-      h1i = 0.5 * (ImagOut[i] - ImagOut[i3]);
-      h2r = 0.5 * (ImagOut[i] + ImagOut[i3]);
-      h2i = -0.5 * (RealOut[i] - RealOut[i3]);
-
-      RealOut[i] = h1r + wr * h2r - wi * h2i;
-      ImagOut[i] = h1i + wr * h2i + wi * h2r;
-      RealOut[i3] = h1r - wr * h2r + wi * h2i;
-      ImagOut[i3] = -h1i + wr * h2i + wi * h2r;
-
-      wr = (wtemp = wr) * wpr - wi * wpi + wr;
-      wi = wi * wpr + wtemp * wpi + wi;
-   }
-
-   RealOut[0] = (h1r = RealOut[0]) + ImagOut[0];
-   ImagOut[0] = h1r - ImagOut[0];
-
-   delete[]tmpReal;
-   delete[]tmpImag;
-#endif //EXPERIMENTAL_USE_REALFFTF
 }
 
-#ifdef EXPERIMENTAL_USE_REALFFTF
 /*
  * InverseRealFFT
  *
@@ -344,10 +274,11 @@ void RealFFT(int NumSamples, float *RealIn, float *RealOut, float *ImagOut)
  * and as a result the output is purely real.
  * Only the first half of RealIn and ImagIn are used due to this
  * symmetry assumption.
+ *
+ * This is merely a wrapper of InverseRealFFTf() from RealFFTf.h.
  */
 void InverseRealFFT(int NumSamples, float *RealIn, float *ImagIn, float *RealOut)
 {
-   // Remap to RealFFTf() function
    int i;
    HFFT hFFT = GetFFT(NumSamples);
    float *pFFT = new float[NumSamples];
@@ -373,24 +304,20 @@ void InverseRealFFT(int NumSamples, float *RealIn, float *ImagIn, float *RealOut
    delete [] pFFT;
    ReleaseFFT(hFFT);
 }
-#endif // EXPERIMENTAL_USE_REALFFTF
 
 /*
  * PowerSpectrum
  *
- * This function computes the same as RealFFT, above, but
- * adds the squares of the real and imaginary part of each
- * coefficient, extracting the power and throwing away the
- * phase.
+ * This function uses RealFFTf() from RealFFTf.h to perform the real
+ * FFT computation, and then squares the real and imaginary part of
+ * each coefficient, extracting the power and throwing away the phase.
  *
  * For speed, it does not call RealFFT, but duplicates some
  * of its code.
  */
 
 void PowerSpectrum(int NumSamples, float *In, float *Out)
 {
-#ifdef EXPERIMENTAL_USE_REALFFTF
-   // Remap to RealFFTf() function
    int i;
    HFFT hFFT = GetFFT(NumSamples);
    float *pFFT = new float[NumSamples];
@@ -411,73 +338,6 @@ void PowerSpectrum(int NumSamples, float *In, float *Out)
    Out[i] = pFFT[1]*pFFT[1];
    delete [] pFFT;
    ReleaseFFT(hFFT);
-
-#else // EXPERIMENTAL_USE_REALFFTF
-
-   int Half = NumSamples / 2;
-   int i;
-
-   float theta = M_PI / Half;
-
-   float *tmpReal = new float[Half];
-   float *tmpImag = new float[Half];
-   float *RealOut = new float[Half];
-   float *ImagOut = new float[Half];
-
-   for (i = 0; i < Half; i++) {
-      tmpReal[i] = In[2 * i];
-      tmpImag[i] = In[2 * i + 1];
-   }
-
-   FFT(Half, 0, tmpReal, tmpImag, RealOut, ImagOut);
-
-   float wtemp = float (sin(0.5 * theta));
-
-   float wpr = -2.0 * wtemp * wtemp;
-   float wpi = -1.0 * float (sin(theta));
-   float wr = 1.0 + wpr;
-   float wi = wpi;
-
-   int i3;
-
-   float h1r, h1i, h2r, h2i, rt, it;
-
-   for (i = 1; i < Half / 2; i++) {
-
-      i3 = Half - i;
-
-      h1r = 0.5 * (RealOut[i] + RealOut[i3]);
-      h1i = 0.5 * (ImagOut[i] - ImagOut[i3]);
-      h2r = 0.5 * (ImagOut[i] + ImagOut[i3]);
-      h2i = -0.5 * (RealOut[i] - RealOut[i3]);
-
-      rt = h1r + wr * h2r - wi * h2i;
-      it = h1i + wr * h2i + wi * h2r;
-
-      Out[i] = rt * rt + it * it;
-
-      rt = h1r - wr * h2r + wi * h2i;
-      it = -h1i + wr * h2i + wi * h2r;
-
-      Out[i3] = rt * rt + it * it;
-
-      wr = (wtemp = wr) * wpr - wi * wpi + wr;
-      wi = wi * wpr + wtemp * wpi + wi;
-   }
-
-   rt = (h1r = RealOut[0]) + ImagOut[0];
-   it = h1r - ImagOut[0];
-   Out[0] = rt * rt + it * it;
-
-   rt = RealOut[Half / 2];
-   it = ImagOut[Half / 2];
-   Out[Half / 2] = rt * rt + it * it;
-
-   delete[]tmpReal;
-   delete[]tmpImag;
-   delete[]RealOut;
-   delete[]ImagOut;
-#endif // EXPERIMENTAL_USE_REALFFTF
 }
 
 /*

src/FFT.h

@@ -77,13 +77,9 @@ void RealFFT(int NumSamples,
  * Computes an Inverse FFT when the input data is conjugate symmetric
  * so the output is purely real.  NumSamples must be a power of
  * two.
- * Requires: EXPERIMENTAL_USE_REALFFTF
  */
-#include "Experimental.h"
-#ifdef EXPERIMENTAL_USE_REALFFTF
 void InverseRealFFT(int NumSamples,
              float *RealIn, float *ImagIn, float *RealOut);
-#endif
 
 /*
  * Computes a FFT of complex input and returns complex output.

src/FreqWindow.cpp

@@ -1255,11 +1255,7 @@ bool SpectrumAnalyst::Calculate(Algorithm alg, int windowFunc,
          case EnhancedAutocorrelation:
 
             // Take FFT
-#ifdef EXPERIMENTAL_USE_REALFFTF
             RealFFT(mWindowSize, in, out, out2);
-#else
-            FFT(mWindowSize, false, in, NULL, out, out2);
-#endif
             // Compute power
             for (int i = 0; i < mWindowSize; i++)
                in[i] = (out[i] * out[i]) + (out2[i] * out2[i]);
@@ -1277,24 +1273,15 @@ bool SpectrumAnalyst::Calculate(Algorithm alg, int windowFunc,
                   in[i] = pow(in[i], 1.0f / 3.0f);
             }
             // Take FFT
-#ifdef EXPERIMENTAL_USE_REALFFTF
             RealFFT(mWindowSize, in, out, out2);
-#else
-            FFT(mWindowSize, false, in, NULL, out, out2);
-#endif
 
             // Take real part of result
             for (int i = 0; i < half; i++)
                mProcessed[i] += out[i];
             break;
 
          case Cepstrum:
-#ifdef EXPERIMENTAL_USE_REALFFTF
             RealFFT(mWindowSize, in, out, out2);
-#else
-            FFT(mWindowSize, false, in, NULL, out, out2);
-#endif
-
             // Compute log power
             // Set a sane lower limit assuming maximum time amplitude of 1.0
             {
@@ -1309,11 +1296,7 @@ bool SpectrumAnalyst::Calculate(Algorithm alg, int windowFunc,
                      in[i] = log(power);
                }
                // Take IFFT
-#ifdef EXPERIMENTAL_USE_REALFFTF
                InverseRealFFT(mWindowSize, in, NULL, out);
-#else
-               FFT(mWindowSize, true, in, NULL, out, out2);
-#endif
 
                // Take real part of result
                for (int i = 0; i < half; i++)

src/Spectrum.cpp

@@ -52,11 +52,7 @@ bool ComputeSpectrum(const float * data, int width,
 
       if (autocorrelation) {
          // Take FFT
-#ifdef EXPERIMENTAL_USE_REALFFTF
          RealFFT(windowSize, in, out, out2);
-#else
-         FFT(windowSize, false, in, NULL, out, out2);
-#endif
          // Compute power
          for (i = 0; i < windowSize; i++)
             in[i] = (out[i] * out[i]) + (out2[i] * out2[i]);
@@ -68,12 +64,7 @@ bool ComputeSpectrum(const float * data, int width,
             in[i] = powf(in[i], 1.0f / 3.0f);
 
          // Take FFT
-#ifdef EXPERIMENTAL_USE_REALFFTF
          RealFFT(windowSize, in, out, out2);
-#else
-         FFT(windowSize, false, in, NULL, out, out2);
-#endif
-
       }
       else
          PowerSpectrum(windowSize, in, out);

src/WaveClip.cpp

@@ -34,6 +34,7 @@
 #include "Resample.h"
 #include "Project.h"
 #include "WaveTrack.h"
+#include "FFT.h"
 
 #include "prefs/SpectrogramSettings.h"
 
@@ -259,8 +260,6 @@ class WaveCache {
 
 };
 
-#ifdef EXPERIMENTAL_USE_REALFFTF
-#include "FFT.h"
 static void ComputeSpectrumUsingRealFFTf
    (float *buffer, HFFT hFFT, const float *window, int len, float *out)
 {
@@ -288,7 +287,6 @@ static void ComputeSpectrumUsingRealFFTf
          out[i] = 10.0*log10f(power);
    }
 }
-#endif // EXPERIMENTAL_USE_REALFFTF
 
 WaveClip::WaveClip(DirManager *projDirManager, sampleFormat format, int rate)
 {
@@ -852,7 +850,6 @@ bool SpecCache::CalculateOneSpectrum
       if (copy)
          useBuffer = scratch;
 
-#ifdef EXPERIMENTAL_USE_REALFFTF
       if (autocorrelation) {
          float *const results = &freq[half * xx];
          // This function does not mutate useBuffer
@@ -956,12 +953,6 @@ bool SpecCache::CalculateOneSpectrum
                results[ii] += gainFactors[ii];
          }
       }
-#else  // EXPERIMENTAL_USE_REALFFTF
-      // This function does not mutate scratch
-      ComputeSpectrum(scratch, windowSize, windowSize,
-         rate, results,
-         autocorrelation, settings.windowType);
-#endif // EXPERIMENTAL_USE_REALFFTF
    }
    return result;
 }
@@ -972,9 +963,7 @@ void SpecCache::Populate
     sampleCount numSamples,
     double offset, double rate, double pixelsPerSecond)
 {
-#ifdef EXPERIMENTAL_USE_REALFFTF
    settings.CacheWindows();
-#endif
 
    const int &frequencyGain = settings.frequencyGain;
    const int &windowSize = settings.windowSize;

src/WaveClip.h

@@ -20,9 +20,7 @@
 #include "xml/XMLTagHandler.h"
 
 #include "Experimental.h"
-#ifdef EXPERIMENTAL_USE_REALFFTF
 #include "RealFFTf.h"
-#endif
 
 #include <wx/gdicmn.h>
 #include <wx/longlong.h>