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