Window function improvemenets: ...

fix off-by-one inconsistencies in the sum-of-cosines windows.

Implement derivatives of the window functions, needed for reassigned
spectrograms.

src/FFT.cpp

@@ -514,82 +514,323 @@ const wxChar *WindowFuncName(int whichFunction)
    }
 }
 
-void WindowFunc(int whichFunction, int NumSamples, float *in)
+void NewWindowFunc(int whichFunction, int NumSamples, bool extraSample, float *in)
 {
-   int i;
+   if (extraSample)
-   double A;
+      --NumSamples;
 
-   switch( whichFunction )
+   switch (whichFunction) {
-   {
    default:
-      fprintf(stderr,"FFT::WindowFunc - Invalid window function: %d\n",whichFunction);
+      fprintf(stderr, "FFT::WindowFunc - Invalid window function: %d\n", whichFunction);
       break;
    case eWinFuncRectangular:
+      // Multiply all by 1.0f -- do nothing
       break;
+
    case eWinFuncBartlett:
+   {
       // Bartlett (triangular) window
-      for (i = 0; i < NumSamples / 2; i++) {
+      const int nPairs = (NumSamples - 1) / 2; // whether even or odd NumSamples, this is correct
-         in[i] *= (i / (float) (NumSamples / 2));
+      const float denom = NumSamples / 2.0f;
-         in[i + (NumSamples / 2)] *=
+      in[0] = 0.0f;
-             (1.0 - (i / (float) (NumSamples / 2)));
+      for (int ii = 1;
+           ii <= nPairs; // Yes, <=
+           ++ii) {
+         const float value = ii / denom;
+         in[ii] *= value;
+         in[NumSamples - ii] *= value;
       }
+      // When NumSamples is even, in[half] should be multiplied by 1.0, so unchanged
+      // When odd, the value of 1.0 is not reached
+   }
       break;
    case eWinFuncHamming:
+   {
       // Hamming
-      for (i = 0; i < NumSamples; i++)
+      const double multiplier = 2 * M_PI / NumSamples;
-         in[i] *= 0.54 - 0.46 * cos(2 * M_PI * i / (NumSamples - 1));
+      static const double coeff0 = 0.54, coeff1 = -0.46;
+      for (int ii = 0; ii < NumSamples; ++ii)
+         in[ii] *= coeff0 + coeff1 * cos(ii * multiplier);
+   }
       break;
    case eWinFuncHanning:
+   {
       // Hanning
-      for (i = 0; i < NumSamples; i++)
+      const double multiplier = 2 * M_PI / NumSamples;
-         in[i] *= 0.50 - 0.50 * cos(2 * M_PI * i / (NumSamples - 1));
+      static const double coeff0 = 0.5, coeff1 = -0.5;
+      for (int ii = 0; ii < NumSamples; ++ii)
+         in[ii] *= coeff0 + coeff1 * cos(ii * multiplier);
+   }
       break;
    case eWinFuncBlackman:
+   {
       // Blackman
-      for (i = 0; i < NumSamples; i++) {
+      const double multiplier = 2 * M_PI / NumSamples;
-          in[i] *= 0.42 - 0.5 * cos (2 * M_PI * i / (NumSamples - 1)) + 0.08 * cos (4 * M_PI * i / (NumSamples - 1));
+      const double multiplier2 = 2 * multiplier;
-      }
+      static const double coeff0 = 0.42, coeff1 = -0.5, coeff2 = 0.08;
+      for (int ii = 0; ii < NumSamples; ++ii)
+         in[ii] *= coeff0 + coeff1 * cos(ii * multiplier) + coeff2 * cos(ii * multiplier2);
+   }
       break;
    case eWinFuncBlackmanHarris:
+   {
       // Blackman-Harris
-      for (i = 0; i < NumSamples; i++) {
+      const double multiplier = 2 * M_PI / NumSamples;
-          in[i] *= 0.35875 - 0.48829 * cos(2 * M_PI * i /(NumSamples-1)) + 0.14128 * cos(4 * M_PI * i/(NumSamples-1)) - 0.01168 * cos(6 * M_PI * i/(NumSamples-1));
+      const double multiplier2 = 2 * multiplier;
-      }
+      const double multiplier3 = 3 * multiplier;
+      static const double coeff0 = 0.35875, coeff1 = -0.48829, coeff2 = 0.14128, coeff3 = -0.01168;
+      for (int ii = 0; ii < NumSamples; ++ii)
+         in[ii] *= coeff0 + coeff1 * cos(ii * multiplier) + coeff2 * cos(ii * multiplier2) + coeff3 * cos(ii * multiplier3);
+   }
       break;
    case eWinFuncWelch:
+   {
       // Welch
-      for (i = 0; i < NumSamples; i++) {
+      const float N = NumSamples;
-          in[i] *= 4*i/(float)NumSamples*(1-(i/(float)NumSamples));
+      for (int ii = 0; ii < NumSamples; ++ii) {
+         const float iOverN = ii / N;
+         in[ii] *= 4 * iOverN * (1 - iOverN);
       }
+   }
+      break;
+   case eWinFuncGaussian25:
+   {
+      // Gaussian (a=2.5)
+      // Precalculate some values, and simplify the fmla to try and reduce overhead
+      static const double A = -2 * 2.5*2.5;
+      const float N = NumSamples;
+      for (int ii = 0; ii < NumSamples; ++ii) {
+         const float iOverN = ii / N;
+         // full
+         // in[ii] *= exp(-0.5*(A*((ii-NumSamples/2)/NumSamples/2))*(A*((ii-NumSamples/2)/NumSamples/2)));
+         // reduced
+         in[ii] *= exp(A * (0.25 + (iOverN * iOverN) - iOverN));
+      }
+   }
+      break;
+   case eWinFuncGaussian35:
+   {
+      // Gaussian (a=3.5)
+      static const double A = -2 * 3.5*3.5;
+      const float N = NumSamples;
+      for (int ii = 0; ii < NumSamples; ++ii) {
+         const float iOverN = ii / N;
+         in[ii] *= exp(A * (0.25 + (iOverN * iOverN) - iOverN));
+      }
+   }
+      break;
+   case eWinFuncGaussian45:
+   {
+      // Gaussian (a=4.5)
+      static const double A = -2 * 4.5*4.5;
+      const float N = NumSamples;
+      for (int ii = 0; ii < NumSamples; ++ii) {
+         const float iOverN = ii / N;
+         in[ii] *= exp(A * (0.25 + (iOverN * iOverN) - iOverN));
+      }
+   }
+      break;
+   }
+
+   if (extraSample && whichFunction != eWinFuncRectangular) {
+      double value = 0.0;
+      switch (whichFunction) {
+      case eWinFuncHamming:
+         value = 0.08;
+         break;
+      case eWinFuncGaussian25:
+         value = exp(-2 * 2.5 * 2.5 * 0.25);
+         break;
+      case eWinFuncGaussian35:
+         value = exp(-2 * 3.5 * 3.5 * 0.25);
+         break;
+      case eWinFuncGaussian45:
+         value = exp(-2 * 4.5 * 4.5 * 0.25);
+         break;
+      default:
+         break;
+      }
+      in[NumSamples] *= value;
+   }
+}
+
+// See cautions in FFT.h !
+void WindowFunc(int whichFunction, int NumSamples, float *in)
+{
+   bool extraSample = false;
+   switch (whichFunction)
+   {
+   case eWinFuncHamming:
+   case eWinFuncHanning:
+   case eWinFuncBlackman:
+   case eWinFuncBlackmanHarris:
+      extraSample = true;
+      break;
+   default:
+      break;
+   case eWinFuncBartlett:
+      // PRL:  Do nothing here either
+      // But I want to comment that the old function did this case
+      // wrong in the second half of the array, in case NumSamples was odd
+      // but I think that never happened, so I am not bothering to preserve that
+      break;
+   }
+   NewWindowFunc(whichFunction, NumSamples, extraSample, in);
+}
+
+void DerivativeOfWindowFunc(int whichFunction, int NumSamples, bool extraSample, float *in)
+{
+   if (eWinFuncRectangular == whichFunction)
+   {
+      // Rectangular
+      // There are deltas at the ends
+      --NumSamples;
+      // in[0] *= 1.0f;
+      for (int ii = 1; ii < NumSamples; ++ii)
+         in[ii] = 0.0f;
+      in[NumSamples] *= -1.0f;
+      return;
+   }
+
+   if (extraSample)
+      --NumSamples;
+
+   double A;
+   switch (whichFunction) {
+   case eWinFuncBartlett:
+   {
+      // Bartlett (triangular) window
+      // There are discontinuities in the derivative at the ends, and maybe at the midpoint
+      const int nPairs = (NumSamples - 1) / 2; // whether even or odd NumSamples, this is correct
+      const float value = 2.0f / NumSamples;
+      in[0] *=
+         // Average the two limiting values of discontinuous derivative
+         value / 2.0f;
+      for (int ii = 1;
+         ii <= nPairs; // Yes, <=
+         ++ii) {
+         in[ii] *= value;
+         in[NumSamples - ii] *= -value;
+      }
+      if (NumSamples % 2 == 0)
+         // Average the two limiting values of discontinuous derivative
+         in[NumSamples / 2] = 0.0f;
+      if (extraSample)
+         in[NumSamples] *=
+         // Average the two limiting values of discontinuous derivative
+            -value / 2.0f;
+      else
+         // Halve the multiplier previously applied
+         // Average the two limiting values of discontinuous derivative
+         in[NumSamples - 1] *= 0.5f;
+   }
+      break;
+   case eWinFuncHamming:
+   {
+      // Hamming
+      // There are deltas at the ends
+      const double multiplier = 2 * M_PI / NumSamples;
+      static const double coeff0 = 0.54, coeff1 = -0.46 * multiplier;
+      in[0] *= coeff0;
+      if (!extraSample)
+         --NumSamples;
+      for (int ii = 0; ii < NumSamples; ++ii)
+         in[ii] *= - coeff1 * sin(ii * multiplier);
+      if (extraSample)
+         in[NumSamples] *= - coeff0;
+      else
+         // slightly different
+         in[NumSamples] *= - coeff0 - coeff1 * sin(NumSamples * multiplier);
+   }
+      break;
+   case eWinFuncHanning:
+   {
+      // Hanning
+      const double multiplier = 2 * M_PI / NumSamples;
+      const double coeff1 = -0.5 * multiplier;
+      for (int ii = 0; ii < NumSamples; ++ii)
+         in[ii] *= - coeff1 * sin(ii * multiplier);
+      if (extraSample)
+         in[NumSamples] = 0.0f;
+   }
+      break;
+   case eWinFuncBlackman:
+   {
+      // Blackman
+      const double multiplier = 2 * M_PI / NumSamples;
+      const double multiplier2 = 2 * multiplier;
+      const double coeff1 = -0.5 * multiplier, coeff2 = 0.08 * multiplier2;
+      for (int ii = 0; ii < NumSamples; ++ii)
+         in[ii] *= - coeff1 * sin(ii * multiplier) - coeff2 * sin(ii * multiplier2);
+      if (extraSample)
+         in[NumSamples] = 0.0f;
+   }
+      break;
+   case eWinFuncBlackmanHarris:
+   {
+      // Blackman-Harris
+      const double multiplier = 2 * M_PI / NumSamples;
+      const double multiplier2 = 2 * multiplier;
+      const double multiplier3 = 3 * multiplier;
+      const double coeff1 = -0.48829 * multiplier,
+         coeff2 = 0.14128 * multiplier2, coeff3 = -0.01168 * multiplier3;
+      for (int ii = 0; ii < NumSamples; ++ii)
+         in[ii] *= - coeff1 * sin(ii * multiplier) - coeff2 * sin(ii * multiplier2) - coeff3 * sin(ii * multiplier3);
+      if (extraSample)
+         in[NumSamples] = 0.0f;
+   }
+      break;
+   case eWinFuncWelch:
+   {
+      // Welch
+      const float N = NumSamples;
+      const float NN = NumSamples * NumSamples;
+      for (int ii = 0; ii < NumSamples; ++ii) {
+         in[ii] *= 4 * (N - ii - ii) / NN;
+      }
+      if (extraSample)
+         in[NumSamples] = 0.0f;
+      // Average the two limiting values of discontinuous derivative
+      in[0] /= 2.0f;
+      in[NumSamples - 1] /= 2.0f;
+   }
       break;
    case eWinFuncGaussian25:
       // Gaussian (a=2.5)
-      // Precalculate some values, and simplify the fmla to try and reduce overhead
+      A = -2 * 2.5*2.5;
-      A=-2*2.5*2.5;
+      goto Gaussian;
-
-      for (i = 0; i < NumSamples; i++) {
-          // full
-          // in[i] *= exp(-0.5*(A*((i-NumSamples/2)/NumSamples/2))*(A*((i-NumSamples/2)/NumSamples/2)));
-          // reduced
-          in[i] *= exp(A*(0.25 + ((i/(float)NumSamples)*(i/(float)NumSamples)) - (i/(float)NumSamples)));
-      }
-      break;
    case eWinFuncGaussian35:
       // Gaussian (a=3.5)
-      A=-2*3.5*3.5;
+      A = -2 * 3.5*3.5;
-      for (i = 0; i < NumSamples; i++) {
+      goto Gaussian;
-          // reduced
-          in[i] *= exp(A*(0.25 + ((i/(float)NumSamples)*(i/(float)NumSamples)) - (i/(float)NumSamples)));
-      }
-      break;
    case eWinFuncGaussian45:
       // Gaussian (a=4.5)
-      A=-2*4.5*4.5;
+      A = -2 * 4.5*4.5;
-
+      goto Gaussian;
-      for (i = 0; i < NumSamples; i++) {
+   Gaussian:
-          // reduced
+   {
-          in[i] *= exp(A*(0.25 + ((i/(float)NumSamples)*(i/(float)NumSamples)) - (i/(float)NumSamples)));
+      // Gaussian (a=2.5)
+      // There are deltas at the ends
+      const float invN = 1.0f / NumSamples;
+      const float invNN = invN * invN;
+      // Simplify formula from the loop for ii == 0, add term for the delta
+      in[0] *= exp(A * 0.25) * (1 - invN);
+      if (!extraSample)
+         --NumSamples;
+      for (int ii = 1; ii < NumSamples; ++ii) {
+         const float iOverN = ii * invN;
+         in[ii] *= exp(A * (0.25 + (iOverN * iOverN) - iOverN)) * (2 * ii * invNN - invN);
       }
+      if (extraSample)
+         in[NumSamples] *= exp(A * 0.25) * (invN - 1);
+      else {
+         // Slightly different
+         const float iOverN = NumSamples * invN;
+         in[NumSamples] *= exp(A * (0.25 + (iOverN * iOverN) - iOverN)) * (2 * NumSamples * invNN - invN - 1);
+      }
+   }
       break;
+   default:
+      fprintf(stderr, "FFT::DerivativeOfWindowFunc - Invalid window function: %d\n", whichFunction);
    }
 }

src/FFT.h

@@ -45,6 +45,9 @@
     * 9: Gaussian(a=4.5)
 */
 
+#include <wx/defs.h>
+#include <wx/wxchar.h>
+
 #ifndef M_PI
 #define	M_PI		3.14159265358979323846  /* pi */
 #endif
@@ -93,18 +96,12 @@ void FFT(int NumSamples,
          float *RealIn, float *ImagIn, float *RealOut, float *ImagOut);
 
 /*
- * Applies a windowing function to the data in place
+ * Multiply values in data by values of the chosen function
- *
+ * DO NOT REUSE!  Prefer NewWindowFunc instead
- * 0: Rectangular (no window)
+ * This version was inconsistent whether the window functions were
- * 1: Bartlett    (triangular)
+ * symmetrical about NumSamples / 2, or about (NumSamples - 1) / 2
- * 2: Hamming
+ * It remains for compatibility until we decide to upgrade all the old uses
- * 3: Hanning
+ * All functions have 0 in data[0] except Rectangular, Hamming and Gaussians
- * 4: Blackman
- * 5: Blackman-Harris
- * 6: Welch
- * 7: Gaussian(a=2.5)
- * 8: Gaussian(a=3.5)
- * 9: Gaussian(a=4.5)
  */
 
 enum eWindowFunctions
@@ -124,6 +121,23 @@ enum eWindowFunctions
 
 void WindowFunc(int whichFunction, int NumSamples, float *data);
 
+/*
+ * Multiply values in data by values of the chosen function
+ * All functions are symmetrical about NumSamples / 2 if extraSample is false,
+ * otherwise about (NumSamples - 1) / 2
+ * All functions have 0 in data[0] except Rectangular, Hamming and Gaussians
+ */
+void NewWindowFunc(int whichFunction, int NumSamples, bool extraSample, float *data);
+
+/*
+ * Multiply values in data by derivative of the chosen function, assuming
+ * sampling interval is unit
+ * All functions are symmetrical about NumSamples / 2 if extraSample is false,
+ * otherwise about (NumSamples - 1) / 2
+ * All functions have 0 in data[0] except Rectangular, Hamming and Gaussians
+ */
+void DerivativeOfWindowFunc(int whichFunction, int NumSamples, bool extraSample, float *data);
+
 /*
  * Returns the name of the windowing function (for UI display)
  */