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:


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.


// Compile with: g++ -std=c++0x -O3 -ffloat-store -fno-unsafe-math-optimizations
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;
  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
        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);
    if (fabs(log(magGoertzel[numDecimated - 1])) >
fabs(log(magMartin[numDecimated - 1]))) {
    if (fabs(log(magGoertzel[numDecimated - 1])) <
fabs(log(magMartin[numDecimated - 1]))) {
  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",
    printf("Martin's RMSE = %.20Lf\n",
  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

Reply via email to