Reassigned spectrograms!! A third algorithm choice in preferences.

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)
  */

src/TrackArtist.cpp

@@ -2375,7 +2375,9 @@ void TrackArtist::DrawClipSpectrum(WaveTrackCache &waveTrackCache,
          (settings, waveTrackCache,
           0, 0, numPixels,
           clip->GetNumSamples(),
-          tOffset, rate);
+          tOffset, rate,
+          0 //FIXME -- make reassignment work with fisheye
+       );
    }
 
    int correctedX = leftOffset - hiddenLeftOffset;

src/WaveClip.cpp

@@ -784,96 +784,200 @@ bool SpecCache::Matches
       algorithm == settings.algorithm;
 }
 
-void SpecCache::CalculateOneSpectrum
+bool SpecCache::CalculateOneSpectrum
    (const SpectrogramSettings &settings,
     WaveTrackCache &waveTrackCache,
     int xx, sampleCount numSamples,
-    double offset, double rate,
+    double offset, double rate, double pixelsPerSecond,
+    int lowerBoundX, int upperBoundX,
     const std::vector<float> &gainFactors,
     float *scratch)
 {
+   bool result = false;
+   const bool reassignment =
+      (settings.algorithm == SpectrogramSettings::algReassignment);
    const int windowSize = settings.windowSize;
-   sampleCount start = where[xx];
+
+   sampleCount start;
+   if (xx < 0)
+      start = where[0] + xx * (rate / pixelsPerSecond);
+   else if (xx > len)
+      start = where[len] + (xx - len) * (rate / pixelsPerSecond);
+   else
+      start = where[xx];
+
    const bool autocorrelation =
       settings.algorithm == SpectrogramSettings::algPitchEAC;
    const int zeroPaddingFactor = (autocorrelation ? 1 : settings.zeroPaddingFactor);
    const int padding = (windowSize * (zeroPaddingFactor - 1)) / 2;
    const int fftLen = windowSize * zeroPaddingFactor;
    const int half = fftLen / 2;
-   float *const results = &freq[half * xx];
-
-   sampleCount len = windowSize;
 
    if (start <= 0 || start >= numSamples) {
-      // Pixel column is out of bounds of the clip!  Should not happen.
+      if (xx >= 0 && xx < len) {
-      std::fill(results, results + half, 0.0f);
+         // Pixel column is out of bounds of the clip!  Should not happen.
+         float *const results = &freq[half * xx];
+         std::fill(results, results + half, 0.0f);
+      }
    }
    else {
-      bool copy = !autocorrelation || (padding > 0);
+      // We can avoid copying memory when ComputeSpectrum is used below
+      bool copy = !autocorrelation || (padding > 0) || reassignment;
       float *useBuffer = 0;
       float *adj = scratch + padding;
 
-      // Take a window of the track centered at this sample.
+      {
-      start -= windowSize >> 1;
+         sampleCount myLen = windowSize;
-      if (start < 0) {
+         // Take a window of the track centered at this sample.
-         // Near the start of the clip, pad left with zeroes as needed.
+         start -= windowSize >> 1;
-         for (sampleCount ii = start; ii < 0; ++ii)
+         if (start < 0) {
-            *adj++ = 0;
+            // Near the start of the clip, pad left with zeroes as needed.
-         len += start;
+            for (sampleCount ii = start; ii < 0; ++ii)
-         start = 0;
+               *adj++ = 0;
-         copy = true;
+            myLen += start;
-      }
+            start = 0;
-      if (start + len > numSamples) {
+            copy = true;
-         // Near the end of the clip, pad right with zeroes as needed.
+         }
-         int newlen = numSamples - start;
-         for (sampleCount ii = newlen; ii < (sampleCount)len; ++ii)
-            adj[ii] = 0;
-         len = newlen;
-         copy = true;
-      }
 
-      if (len > 0) {
+         if (start + myLen > numSamples) {
-         // Copy samples out of the track.
+            // Near the end of the clip, pad right with zeroes as needed.
-         useBuffer = (float*)(waveTrackCache.Get(floatSample,
+            int newlen = numSamples - start;
-                             floor(0.5 + start + offset * rate), len));
+            for (sampleCount ii = newlen; ii < (sampleCount)myLen; ++ii)
-         if (copy)
+               adj[ii] = 0;
-            memcpy(adj, useBuffer, len * sizeof(float));
+            myLen = newlen;
+            copy = true;
+         }
+
+         if (myLen > 0) {
+            useBuffer = (float*)(waveTrackCache.Get(floatSample,
+               floor(0.5 + start + offset * rate), myLen));
+            if (copy)
+               memcpy(adj, useBuffer, myLen * sizeof(float));
+         }
       }
 
       if (copy)
          useBuffer = scratch;
 
 #ifdef EXPERIMENTAL_USE_REALFFTF
-      if (autocorrelation)
+      if (autocorrelation) {
+         float *const results = &freq[half * xx];
+         // This function does not mutate useBuffer
          ComputeSpectrum(useBuffer, windowSize, windowSize,
             rate, results,
             autocorrelation, settings.windowType);
-      else
+      }
-         // Do the FFT.  Note that scratch is multiplied by the window,
+      else if (reassignment) {
+         static const double epsilon = 1e-16;
+         const HFFT hFFT = settings.hFFT;
+
+         float *const scratch2 = scratch + fftLen;
+         std::copy(scratch, scratch2, scratch2);
+
+         float *const scratch3 = scratch + 2 * fftLen;
+         std::copy(scratch, scratch2, scratch3);
+
+         {
+            const float *const window = settings.window;
+            for (int ii = 0; ii < fftLen; ++ii)
+               scratch[ii] *= window[ii];
+            RealFFTf(scratch, hFFT);
+         }
+
+         {
+            const float *const dWindow = settings.dWindow;
+            for (int ii = 0; ii < fftLen; ++ii)
+               scratch2[ii] *= dWindow[ii];
+            RealFFTf(scratch2, hFFT);
+         }
+
+         {
+            const float *const tWindow = settings.tWindow;
+            for (int ii = 0; ii < fftLen; ++ii)
+               scratch3[ii] *= tWindow[ii];
+            RealFFTf(scratch3, hFFT);
+         }
+
+         for (int ii = 0; ii < hFFT->Points; ++ii) {
+            const int index = hFFT->BitReversed[ii];
+            const float
+               denomRe = scratch[index],
+               denomIm = ii == 0 ? 0 : scratch[index + 1];
+            const double power = denomRe * denomRe + denomIm * denomIm;
+            if (power < epsilon)
+               // Avoid dividing by near-zero below
+               continue;
+
+            double freqCorrection;
+            {
+               const double multiplier = -fftLen / (2.0f * M_PI);
+               const float
+                  numRe = scratch2[index],
+                  numIm = ii == 0 ? 0 : scratch2[index + 1];
+               // Find complex quotient --
+               // Which means, multiply numerator by conjugate of denominator,
+               // then divide by norm squared of denominator --
+               // Then just take its imaginary part.
+               const double
+                  quotIm = (-numRe * denomIm + numIm * denomRe) / power;
+               // With appropriate multiplier, that becomes the correction of
+               // the frequency bin.
+               freqCorrection = multiplier * quotIm;
+            }
+
+            const int bin = int(ii + freqCorrection + 0.5f);
+            if (bin >= 0 && bin < hFFT->Points) {
+               double timeCorrection;
+               {
+                  const float
+                     numRe = scratch3[index],
+                     numIm = ii == 0 ? 0 : scratch3[index + 1];
+                  // Find another complex quotient --
+                  // Then just take its real part.
+                  // The result has sample interval as unit.
+                  timeCorrection =
+                     (numRe * denomRe + numIm * denomIm) / power;
+               }
+
+               int correctedX = (floor(0.5 + xx + timeCorrection * pixelsPerSecond / rate));
+               if (correctedX >= lowerBoundX && correctedX < upperBoundX)
+                  result = true,
+                  freq[half * correctedX + bin] += power;
+            }
+         }
+      }
+      else {
+         float *const results = &freq[half * xx];
+
+         // Do the FFT.  Note that useBuffer is multiplied by the window,
          // and the window is initialized with leading and trailing zeroes
          // when there is padding.  Therefore we did not need to reinitialize
-         // the part of scratch in the padding zones.
+         // the part of useBuffer in the padding zones.
+
+         // This function mutates useBuffer
          ComputeSpectrumUsingRealFFTf
             (useBuffer, settings.hFFT, settings.window, fftLen, results);
+         if (!gainFactors.empty()) {
+            // Apply a frequency-dependant gain factor
+            for (int ii = 0; ii < half; ++ii)
+               results[ii] += gainFactors[ii];
+         }
+      }
 #else  // EXPERIMENTAL_USE_REALFFTF
-      ComputeSpectrum(buffer, windowSize, windowSize,
+      // This function does not mutate scratch
+      ComputeSpectrum(scratch, windowSize, windowSize,
          rate, results,
          autocorrelation, settings.windowType);
 #endif // EXPERIMENTAL_USE_REALFFTF
-      if (!autocorrelation &&
-          !gainFactors.empty()) {
-         // Apply a frequency-dependant gain factor
-         for (int ii = 0; ii < half; ++ii)
-            results[ii] += gainFactors[ii];
-      }
    }
+   return result;
 }
 
 void SpecCache::Populate
    (const SpectrogramSettings &settings, WaveTrackCache &waveTrackCache,
     int copyBegin, int copyEnd, int numPixels,
     sampleCount numSamples,
-    double offset, double rate)
+    double offset, double rate, double pixelsPerSecond)
 {
 #ifdef EXPERIMENTAL_USE_REALFFTF
    settings.CacheWindows();
@@ -883,6 +987,8 @@ void SpecCache::Populate
    const int &windowSize = settings.windowSize;
    const bool autocorrelation =
       settings.algorithm == SpectrogramSettings::algPitchEAC;
+   const bool reassignment =
+      settings.algorithm == SpectrogramSettings::algReassignment;
 #ifdef EXPERIMENTAL_ZERO_PADDED_SPECTROGRAMS
    const int &zeroPaddingFactor = autocorrelation ? 1 : settings.zeroPaddingFactor;
 #else
@@ -892,13 +998,15 @@ void SpecCache::Populate
    // FFT length may be longer than the window of samples that affect results
    // because of zero padding done for increased frequency resolution
    const int fftLen = windowSize * zeroPaddingFactor;
+   const int half = fftLen / 2;
 
-   std::vector<float> buffer(
+   const size_t bufferSize = fftLen;
-      fftLen
+
-   );
+   std::vector<float> buffer(reassignment ? 3 * bufferSize : bufferSize);
 
    std::vector<float> gainFactors;
-   ComputeSpectrogramGainFactors(fftLen, rate, frequencyGain, gainFactors);
+   if (!autocorrelation)
+      ComputeSpectrogramGainFactors(fftLen, rate, frequencyGain, gainFactors);
 
    // Loop over the ranges before and after the copied portion and compute anew.
    // One of the ranges may be empty.
@@ -907,8 +1015,62 @@ void SpecCache::Populate
       const int upperBoundX = jj == 0 ? copyBegin : numPixels;
       for (sampleCount xx = lowerBoundX; xx < upperBoundX; ++xx)
          CalculateOneSpectrum(
-         settings, waveTrackCache, xx, numSamples,
+            settings, waveTrackCache, xx, numSamples,
-         offset, rate, gainFactors, &buffer[0]);
+            offset, rate, pixelsPerSecond,
+            lowerBoundX, upperBoundX,
+            gainFactors, &buffer[0]);
+
+      if (reassignment) {
+         // Need to look beyond the edges of the range to accumulate more
+         // time reassignments.
+         // I'm not sure what's a good stopping criterion?
+         sampleCount xx = lowerBoundX;
+         const double pixelsPerSample = pixelsPerSecond / rate;
+         const int limit = std::min(int(0.5 + fftLen * pixelsPerSample), 100);
+         for (int ii = 0; ii < limit; ++ii)
+         {
+            const bool result =
+               CalculateOneSpectrum(
+                  settings, waveTrackCache, --xx, numSamples,
+                  offset, rate, pixelsPerSecond,
+                  lowerBoundX, upperBoundX,
+                  gainFactors, &buffer[0]);
+            if (!result)
+               break;
+         }
+
+         xx = upperBoundX;
+         for (int ii = 0; ii < limit; ++ii)
+         {
+            const bool result =
+               CalculateOneSpectrum(
+                  settings, waveTrackCache, xx++, numSamples,
+                  offset, rate, pixelsPerSecond,
+                  lowerBoundX, upperBoundX,
+                  gainFactors, &buffer[0]);
+            if (!result)
+               break;
+         }
+
+         // Now Convert to dB terms.  Do this only after accumulating
+         // power values, which may cross columns with the time correction.
+         for (sampleCount xx = lowerBoundX; xx < upperBoundX; ++xx) {
+            float *const results = &freq[half * xx];
+            const HFFT hFFT = settings.hFFT;
+            for (int ii = 0; ii < hFFT->Points; ++ii) {
+               float &power = results[ii];
+               if (power <= 0)
+                  power = -160.0;
+               else
+                  power = 10.0*log10f(power);
+            }
+            if (!gainFactors.empty()) {
+               // Apply a frequency-dependant gain factor
+               for (int ii = 0; ii < half; ++ii)
+                  results[ii] += gainFactors[ii];
+            }
+         }
+      }
    }
 }
 
@@ -935,7 +1097,7 @@ bool WaveClip::GetSpectrogram(WaveTrackCache &waveTrackCache,
    const int fftLen = windowSize * zeroPaddingFactor;
    const int half = fftLen / 2;
 
-   const bool match =
+   bool match =
       mSpecCache &&
       mSpecCache->len > 0 &&
       mSpecCache->Matches
@@ -949,6 +1111,11 @@ bool WaveClip::GetSpectrogram(WaveTrackCache &waveTrackCache,
       return false;  //hit cache completely
    }
 
+   if (settings.algorithm == SpectrogramSettings::algReassignment)
+      // Caching is not implemented for reassignment, unless for
+      // a complete hit, because of the complications of time reassignment
+      match = false;
+
    std::auto_ptr<SpecCache> oldCache(mSpecCache);
    mSpecCache = 0;
 
@@ -979,6 +1146,8 @@ bool WaveClip::GetSpectrogram(WaveTrackCache &waveTrackCache,
       numPixels, settings.algorithm, pixelsPerSecond, t0,
       windowType, windowSize, zeroPaddingFactor, frequencyGain);
 
+   // purposely offset the display 1/2 sample to the left (as compared
+   // to waveform display) to properly center response of the FFT
    fillWhere(mSpecCache->where, numPixels, 0.5, correction,
       t0, mRate, samplesPerPixel);
 
@@ -993,7 +1162,8 @@ bool WaveClip::GetSpectrogram(WaveTrackCache &waveTrackCache,
 
    mSpecCache->Populate
       (settings, waveTrackCache, copyBegin, copyEnd, numPixels,
-       mSequence->GetNumSamples(), mOffset, mRate);
+       mSequence->GetNumSamples(),
+       mOffset, mRate, pixelsPerSecond);
 
    mSpecCache->dirty = mDirty;
    spectrogram = &mSpecCache->freq[0];

src/WaveClip.h

@@ -92,11 +92,12 @@ public:
    bool Matches(int dirty_, double pixelsPerSecond,
       const SpectrogramSettings &settings, double rate) const;
 
-   void CalculateOneSpectrum
+   bool CalculateOneSpectrum
       (const SpectrogramSettings &settings,
        WaveTrackCache &waveTrackCache,
        int xx, sampleCount numSamples,
-       double offset, double rate,
+       double offset, double rate, double pixelsPerSecond,
+       int lowerBoundX, int upperBoundX,
        const std::vector<float> &gainFactors,
        float *scratch);
 
@@ -104,7 +105,7 @@ public:
       (const SpectrogramSettings &settings, WaveTrackCache &waveTrackCache,
        int copyBegin, int copyEnd, int numPixels,
        sampleCount numSamples,
-       double offset, double rate);
+       double offset, double rate, double pixelsPerSecond);
 
    const int          len; // counts pixels, not samples
    const int          algorithm;

src/prefs/SpectrogramSettings.cpp

@@ -56,6 +56,8 @@ SpectrogramSettings::Globals
 SpectrogramSettings::SpectrogramSettings()
    : hFFT(0)
    , window(0)
+   , dWindow(0)
+   , tWindow(0)
 {
    LoadPrefs();
 }
@@ -92,6 +94,8 @@ SpectrogramSettings::SpectrogramSettings(const SpectrogramSettings &other)
    // Do not copy these!
    , hFFT(0)
    , window(0)
+   , tWindow(0)
+   , dWindow(0)
 {
 }
 
@@ -190,6 +194,8 @@ const wxArrayString &SpectrogramSettings::GetAlgorithmNames()
    if (theArray.IsEmpty()) {
       // Keep in correspondence with enum SpectrogramSettings::Algorithm:
       theArray.Add(_("Frequencies"));
+      /* i18n-hint: the Reassignment algorithm for spectrograms */
+      theArray.Add(_("Reassignment"));
       /* i18n-hint: EAC abbreviates "Enhanced Autocorrelation" */
       theArray.Add(_("Pitch (EAC)"));
    }
@@ -386,6 +392,14 @@ void SpectrogramSettings::DestroyWindows()
       delete[] window;
       window = NULL;
    }
+   if (dWindow != NULL) {
+      delete[] dWindow;
+      dWindow = NULL;
+   }
+   if (tWindow != NULL) {
+      delete[] tWindow;
+      tWindow = NULL;
+   }
 #endif
 }
 
@@ -403,7 +417,11 @@ namespace
       window = new float[fftLen];
       int ii;
 
+      const bool extra = padding > 0;
       wxASSERT(windowSize % 2 == 0);
+      if (extra)
+         // For windows that do not go to 0 at the edges, this improves symmetry
+         ++windowSize;
       const int endOfWindow = padding + windowSize;
       // Left and right padding
       for (ii = 0; ii < padding; ++ii) {
@@ -416,25 +434,17 @@ namespace
       // Overwrite middle as needed
       switch (which) {
       case WINDOW:
-         WindowFunc(windowType, windowSize, window + padding);
+         NewWindowFunc(windowType, windowSize, extra, window + padding);
-         // NewWindowFunc(windowType, windowSize, extra, window + padding);
          break;
-      case TWINDOW:
-         wxASSERT(false);
-#if 0
          // Future, reassignment
+      case TWINDOW:
          NewWindowFunc(windowType, windowSize, extra, window + padding);
          for (int ii = padding, multiplier = -windowSize / 2; ii < endOfWindow; ++ii, ++multiplier)
             window[ii] *= multiplier;
          break;
-#endif
       case DWINDOW:
-         wxASSERT(false);
-#if 0
-         // Future, reassignment
          DerivativeOfWindowFunc(windowType, windowSize, extra, window + padding);
          break;
-#endif
       default:
          wxASSERT(false);
       }
@@ -464,6 +474,10 @@ void SpectrogramSettings::CacheWindows() const
          EndFFT(hFFT);
       hFFT = InitializeFFT(fftLen);
       RecreateWindow(window, WINDOW, fftLen, padding, windowType, windowSize, scale);
+      if (algorithm == algReassignment) {
+         RecreateWindow(tWindow, TWINDOW, fftLen, padding, windowType, windowSize, scale);
+         RecreateWindow(dWindow, DWINDOW, fftLen, padding, windowType, windowSize, scale);
+      }
    }
 #endif // EXPERIMENTAL_USE_REALFFTF
 }
@@ -557,7 +571,7 @@ int SpectrogramSettings::GetFFTLength() const
 {
    return windowSize
 #ifdef EXPERIMENTAL_ZERO_PADDED_SPECTROGRAMS
-      * ((algorithm == algSTFT) ? zeroPaddingFactor : 1);
+      * ((algorithm != algPitchEAC) ? zeroPaddingFactor : 1);
 #endif
    ;
 }

src/prefs/SpectrogramSettings.h

@@ -130,6 +130,7 @@ public:
 
    enum Algorithm {
       algSTFT = 0,
+      algReassignment,
       algPitchEAC,
 
       algNumAlgorithms,
@@ -153,6 +154,11 @@ public:
    // Variables used for computing the spectrum
    mutable FFTParam      *hFFT;
    mutable float         *window;
+
+   // Two other windows for computing reassigned spectrogram
+   mutable float         *tWindow; // Window times time parameter
+   mutable float         *dWindow; // Derivative of window
+
 #endif
 };
 #endif

src/prefs/SpectrumPrefs.cpp

@@ -459,7 +459,8 @@ void SpectrumPrefs::OnAlgorithm(wxCommandEvent &evt)
 void SpectrumPrefs::EnableDisableSTFTOnlyControls()
 {
    // Enable or disable other controls that are applicable only to STFT.
-   const bool STFT = (mAlgorithmChoice->GetSelection() == 0);
+   const bool STFT =
+      (mAlgorithmChoice->GetSelection() != SpectrogramSettings::algPitchEAC);
    mGain->Enable(STFT);
    mRange->Enable(STFT);
    mFrequencyGain->Enable(STFT);