On Thu, Feb 21, 2019 at 11:16 PM robert bristow-johnson <r...@audioimagination.com> wrote:
> But Martin, if you let this thing run for days on end, would not eventually > the amplitude of the output change a bit? Short answer: yes, sometimes significantly for audio purposes when using 32-bit float state variables, but not significantly when using 64-bit doubles. Even as poles lie exactly on the unit circle, limited precision in the state variables leads to rounding errors. Because of the recursive nature of the algorithm, rounding errors accumulate in manner of a (pseudo)random walk. Almost always a past state will be encountered again, and the system enters a limit cycle, locking the error into a range of error values close to what was accumulated thus far. With 32-bit floating point state variables, sometimes significant numerical error is accumulated before a limit cycle is entered, depending on pseudo-luck with the sinusoid frequency. With 64-bit floating point state variables the problem seems to be almost never significant even though entering a limit cycle becomes a rare occasion. The error accumulation and limit cycles can be visualized by plotting the magnitude of the quadrature pair output of Martin's algorithm (red), and the same for a Goertzel algorithm based quadrature oscillator (blue) for comparison, as function of output sample number: https://imgur.com/a/cE5R50d For this plot, 32-bit float variables were used, and rounding mode was set to round-to-nearest, and angular frequency arbitrarily to (float) sqrt(2.0) radians to get visible and somewhat lengthy limit cycles. For the Goertzel algorithm, about 4 periods of a limit cycle, and for Martin's algorithm, about 2 periods of a limit cycle are visible, for the first one million output samples. Based on testing with random frequencies and single (32-bit) and double (64-bit) precision floating point numbers, compared to Goertzel, Martin's algorithm tends to have faster error accumulation (due to a longer cascade of rounded operations per iteration), giving a larger error, measured as absolute log magnitude, for about 62 % of random frequencies, both for floats and doubles (tested over 10^4 replications, measured at 10^4 samples generated). At one million samples generated, the root mean square magnitude error was roughly 0.0002 for Goertzel and 0.002 for Martin's algorithm using single precision state variables, and roughly 7*10^-13 for Goertzel and 3*10^-12 for Martin's algorithm using double precision floating point state variables (tested over 1000 replications). C++ test code below. -olli // Compile with: g++ -std=c++0x -O3 -ffloat-store -fno-unsafe-math-optimizations #define _USE_MATH_DEFINES typedef float FLOAT_TYPE; // float or double const bool printStats = true; // Print statistics const bool printData = false; // Print (last) data #include <math.h> #include <stdio.h> #include <fenv.h> #include <random> #include <chrono> main() { long int numReps = 1000; long int numSamples = 1000; long int decimate = 1; // Decimate results by this ratio (useful if huge numReps) unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); std::mt19937 rng(seed); std::uniform_real_distribution<long double> uniform(0.0, 1.0); int numDecimated = numSamples/decimate + 1; long double *magGoertzel = new long double[numDecimated]; long double *magMartin = new long double[numDecimated]; long double *sseGoertzel = new long double[numDecimated]; long double *sseMartin = new long double[numDecimated]; for (long int count = 0; count < numDecimated; count++) { sseGoertzel[count] = 0; sseMartin[count] = 0; } fesetround(FE_TONEAREST); // FE_TONEAREST, FE_UPWARD, FE_DOWNWARD, or FE_TOWARDZERO long int numGoertzelWorse = 0; long int numMartinWorse = 0; for (long int rep = 0; rep < numReps; rep++) { FLOAT_TYPE w = M_PI*uniform(rng); // Goertzel FLOAT_TYPE coef = 2*cos(w); FLOAT_TYPE outHistory1 = 0, outHistory2 = 0; FLOAT_TYPE input = 1; // First input value FLOAT_TYPE cosw = cos(w); FLOAT_TYPE sinw = sin(w); // Martin FLOAT_TYPE k1 = tan(w/2); FLOAT_TYPE k2 = 2*k1/(1 + k1*k1); FLOAT_TYPE u = 1; FLOAT_TYPE v = 0; numDecimated = 0; for (long int count = 0; count < numSamples; count++) { // Goertzel FLOAT_TYPE out = input + coef*outHistory1 - outHistory2; outHistory2 = outHistory1; outHistory1 = out; input = 0; FLOAT_TYPE re = outHistory1 - outHistory2*cosw; FLOAT_TYPE im = outHistory2*sinw; // Martin FLOAT_TYPE tmp = u - k1*v; v = v + k2*tmp; u = tmp - k1*v; // Results if (count % decimate == 0) { magGoertzel[numDecimated] = sqrt(re*(long double)re + im*(long double)im); magMartin[numDecimated] = (FLOAT_TYPE)sqrt(u*(long double)u + v*(long double)v); sseGoertzel[numDecimated] += pow(magGoertzel[numDecimated] - 1, 2); sseMartin[numDecimated] += pow(magMartin[numDecimated] - 1, 2); numDecimated++; } } if (fabs(log(magGoertzel[numDecimated - 1])) > fabs(log(magMartin[numDecimated - 1]))) { numGoertzelWorse++; } if (fabs(log(magGoertzel[numDecimated - 1])) < fabs(log(magMartin[numDecimated - 1]))) { numMartinWorse++; } } if (printStats) { printf("Goertzel worse ratio = %f\n", numGoertzelWorse/(FLOAT_TYPE)numReps); printf("Martin's worse ratio = %f\n", numMartinWorse/(FLOAT_TYPE)numReps); printf("Goertzel RMSE = %.20Lf\n", sqrtl(sseGoertzel[numDecimated-1]/numReps)); printf("Martin's RMSE = %.20Lf\n", sqrtl(sseMartin[numDecimated-1]/numReps)); } if (printData) { for (int count = 0; count < numDecimated; count++) { printf("%d,%.20Lf,%.20Lf,%.20Lf,%.20Lf\n", count, magGoertzel[count]-1, magMartin[count]-1, sqrtl(sseGoertzel[count]/numReps), sqrtl(sseMartin[count]/numReps)); } } } _______________________________________________ dupswapdrop: music-dsp mailing list music-dsp@music.columbia.edu https://lists.columbia.edu/mailman/listinfo/music-dsp
