diff --git a/src/FFT.cpp b/src/FFT.cpp index af1217bf7..00922d2bc 100644 --- a/src/FFT.cpp +++ b/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; - double A; + if (extraSample) + --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++) { - in[i] *= (i / (float) (NumSamples / 2)); - in[i + (NumSamples / 2)] *= - (1.0 - (i / (float) (NumSamples / 2))); + const int nPairs = (NumSamples - 1) / 2; // whether even or odd NumSamples, this is correct + const float denom = NumSamples / 2.0f; + in[0] = 0.0f; + 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++) - in[i] *= 0.54 - 0.46 * cos(2 * M_PI * i / (NumSamples - 1)); + const double multiplier = 2 * M_PI / NumSamples; + 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++) - in[i] *= 0.50 - 0.50 * cos(2 * M_PI * i / (NumSamples - 1)); + const double multiplier = 2 * M_PI / NumSamples; + 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++) { - in[i] *= 0.42 - 0.5 * cos (2 * M_PI * i / (NumSamples - 1)) + 0.08 * cos (4 * M_PI * i / (NumSamples - 1)); - } + const double multiplier = 2 * M_PI / NumSamples; + 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++) { - 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 multiplier = 2 * M_PI / NumSamples; + 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++) { - in[i] *= 4*i/(float)NumSamples*(1-(i/(float)NumSamples)); + const float N = 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; - - 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; + A = -2 * 2.5*2.5; + goto Gaussian; case eWinFuncGaussian35: // Gaussian (a=3.5) - A=-2*3.5*3.5; - for (i = 0; i < NumSamples; i++) { - // reduced - in[i] *= exp(A*(0.25 + ((i/(float)NumSamples)*(i/(float)NumSamples)) - (i/(float)NumSamples))); - } - break; + A = -2 * 3.5*3.5; + goto Gaussian; case eWinFuncGaussian45: // Gaussian (a=4.5) - A=-2*4.5*4.5; - - for (i = 0; i < NumSamples; i++) { - // reduced - in[i] *= exp(A*(0.25 + ((i/(float)NumSamples)*(i/(float)NumSamples)) - (i/(float)NumSamples))); + A = -2 * 4.5*4.5; + goto Gaussian; + Gaussian: + { + // 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); } } diff --git a/src/FFT.h b/src/FFT.h index 80de2b398..10dad8053 100644 --- a/src/FFT.h +++ b/src/FFT.h @@ -45,6 +45,9 @@ * 9: Gaussian(a=4.5) */ +#include +#include + #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 - * - * 0: Rectangular (no window) - * 1: Bartlett (triangular) - * 2: Hamming - * 3: Hanning - * 4: Blackman - * 5: Blackman-Harris - * 6: Welch - * 7: Gaussian(a=2.5) - * 8: Gaussian(a=3.5) - * 9: Gaussian(a=4.5) + * Multiply values in data by values of the chosen function + * DO NOT REUSE! Prefer NewWindowFunc instead + * This version was inconsistent whether the window functions were + * symmetrical about NumSamples / 2, or about (NumSamples - 1) / 2 + * It remains for compatibility until we decide to upgrade all the old uses + * All functions have 0 in data[0] except Rectangular, Hamming and Gaussians */ 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) */ diff --git a/src/TrackArtist.cpp b/src/TrackArtist.cpp index 802f066a3..e8fb61653 100644 --- a/src/TrackArtist.cpp +++ b/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; diff --git a/src/WaveClip.cpp b/src/WaveClip.cpp index 4146bf3ed..82a95b409 100644 --- a/src/WaveClip.cpp +++ b/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 &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. - std::fill(results, results + half, 0.0f); + if (xx >= 0 && xx < len) { + // 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; - if (start < 0) { - // Near the start of the clip, pad left with zeroes as needed. - for (sampleCount ii = start; ii < 0; ++ii) - *adj++ = 0; - len += start; - start = 0; - copy = true; - } - if (start + len > numSamples) { - // 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; - } + { + sampleCount myLen = windowSize; + // Take a window of the track centered at this sample. + start -= windowSize >> 1; + if (start < 0) { + // Near the start of the clip, pad left with zeroes as needed. + for (sampleCount ii = start; ii < 0; ++ii) + *adj++ = 0; + myLen += start; + start = 0; + copy = true; + } - if (len > 0) { - // Copy samples out of the track. - useBuffer = (float*)(waveTrackCache.Get(floatSample, - floor(0.5 + start + offset * rate), len)); - if (copy) - memcpy(adj, useBuffer, len * sizeof(float)); + if (start + myLen > numSamples) { + // Near the end of the clip, pad right with zeroes as needed. + int newlen = numSamples - start; + for (sampleCount ii = newlen; ii < (sampleCount)myLen; ++ii) + adj[ii] = 0; + 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 buffer( - fftLen - ); + const size_t bufferSize = fftLen; + + std::vector buffer(reassignment ? 3 * bufferSize : bufferSize); std::vector 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, - offset, rate, gainFactors, &buffer[0]); + settings, waveTrackCache, xx, numSamples, + 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 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]; diff --git a/src/WaveClip.h b/src/WaveClip.h index 09de9c8f0..35dc76fc5 100644 --- a/src/WaveClip.h +++ b/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 &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; diff --git a/src/prefs/SpectrogramSettings.cpp b/src/prefs/SpectrogramSettings.cpp index 4b8ed4985..ef0b6227d 100644 --- a/src/prefs/SpectrogramSettings.cpp +++ b/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 ; } diff --git a/src/prefs/SpectrogramSettings.h b/src/prefs/SpectrogramSettings.h index 57c1cc58a..5077a188f 100644 --- a/src/prefs/SpectrogramSettings.h +++ b/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 diff --git a/src/prefs/SpectrumPrefs.cpp b/src/prefs/SpectrumPrefs.cpp index b662a6f12..5e125b530 100644 --- a/src/prefs/SpectrumPrefs.cpp +++ b/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);