reassignment, internals, implementing time correction, but -- ...

... it does not yet interact correctly with caching, so results may be
slightly wrong if you scroll the view by less than a screen.

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

@@ -788,7 +788,8 @@ void 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)
 {
@@ -864,6 +865,9 @@ void SpecCache::CalculateOneSpectrum
          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)
@@ -878,7 +882,13 @@ void SpecCache::CalculateOneSpectrum
             RealFFTf(scratch2, hFFT);
          }
 
-         const double multiplier = -fftLen / (2.0f * M_PI);
+         {
+            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
@@ -889,30 +899,41 @@ void SpecCache::CalculateOneSpectrum
                // Avoid dividing by near-zero below
                continue;
 
-            const float
+            double freqCorrection;
-               numRe = scratch2[index],
+            {
-               numIm = ii == 0 ? 0 : scratch2[index + 1];
+               const double multiplier = -fftLen / (2.0f * M_PI);
-            // Find complex quotient --
+               const float
-            // Which means, multiply numerator by conjugate of denominator,
+                  numRe = scratch2[index],
-            // then divide by norm squared of denominator --
+                  numIm = ii == 0 ? 0 : scratch2[index + 1];
-            // Then just take its imaginary part.
+               // Find complex quotient --
-            const double
+               // Which means, multiply numerator by conjugate of denominator,
-               quotIm = (-numRe * denomIm + numIm * denomRe) / power;
+               // then divide by norm squared of denominator --
-            // With appropriate multiplier, that becomes the correction of
+               // Then just take its imaginary part.
-            // the frequency bin.
+               const double
-            const double correction = multiplier * quotIm;
+                  quotIm = (-numRe * denomIm + numIm * denomRe) / power;
-            const int bin = int(ii + correction + 0.5f);
+               // With appropriate multiplier, that becomes the correction of
-            if (bin >= 0 && bin < hFFT->Points)
+               // the frequency bin.
-               freq[half * xx + bin] += power;
+               freqCorrection = multiplier * quotIm;
-         }
+            }
 
-         // Now Convert to dB terms
+            const int bin = int(ii + freqCorrection + 0.5f);
-         for (int ii = 0; ii < hFFT->Points; ++ii) {
+            if (bin >= 0 && bin < hFFT->Points) {
-            float &power = freq[half * xx + ii];
+               double timeCorrection;
-            if (power <= 0)
+               {
-               power = -160.0;
+                  const float
-            else
+                     numRe = scratch3[index],
-               power = 10.0*log10f(power);
+                     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)
+                  freq[half * correctedX + bin] += power;
+            }
          }
       }
       else {
@@ -945,7 +966,7 @@ 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();
@@ -970,7 +991,7 @@ void SpecCache::Populate
 
    const size_t bufferSize = fftLen;
 
-   std::vector<float> buffer(reassignment ? 2 * bufferSize : bufferSize);
+   std::vector<float> buffer(reassignment ? 3 * bufferSize : bufferSize);
 
    std::vector<float> gainFactors;
    if (!autocorrelation)
@@ -983,8 +1004,10 @@ 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) {
          // Now Convert to dB terms.  Do this only after accumulating
@@ -1092,7 +1115,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

@@ -96,7 +96,8 @@ public:
       (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

@@ -57,6 +57,7 @@ SpectrogramSettings::SpectrogramSettings()
    : hFFT(0)
    , window(0)
    , dWindow(0)
+   , tWindow(0)
 {
    LoadPrefs();
 }
@@ -93,9 +94,7 @@ SpectrogramSettings::SpectrogramSettings(const SpectrogramSettings &other)
    // Do not copy these!
    , hFFT(0)
    , window(0)
-#if 0
    , tWindow(0)
-#endif
    , dWindow(0)
 {
 }
@@ -397,6 +396,10 @@ void SpectrogramSettings::DestroyWindows()
       delete[] dWindow;
       dWindow = NULL;
    }
+   if (tWindow != NULL) {
+      delete[] tWindow;
+      tWindow = NULL;
+   }
 #endif
 }
 
@@ -433,14 +436,12 @@ namespace
       case WINDOW:
          NewWindowFunc(windowType, windowSize, extra, window + padding);
          break;
-#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:
          DerivativeOfWindowFunc(windowType, windowSize, extra, window + padding);
          break;
@@ -474,9 +475,7 @@ void SpectrogramSettings::CacheWindows() const
       hFFT = InitializeFFT(fftLen);
       RecreateWindow(window, WINDOW, fftLen, padding, windowType, windowSize, scale);
       if (algorithm == algReassignment) {
-#if 0
          RecreateWindow(tWindow, TWINDOW, fftLen, padding, windowType, windowSize, scale);
-#endif
          RecreateWindow(dWindow, DWINDOW, fftLen, padding, windowType, windowSize, scale);
       }
    }

src/prefs/SpectrogramSettings.h

@@ -156,9 +156,7 @@ public:
    mutable float         *window;
 
    // Two other windows for computing reassigned spectrogram
-#if 0
    mutable float         *tWindow; // Window times time parameter
-#endif
    mutable float         *dWindow; // Derivative of window
 
 #endif