(* Works, missing some of the output from the other implementations *)
MODULE Goertzel;
  IMPORT Texts, Oberon, Math;
  CONST SamplingRate* = 8000.0;
    TargetFrequency* = 941.0;
    N* = 205;
    pi* = 3.141592654;

  TYPE Sample* = BYTE;

  VAR coeff, q1, q2, sine, cosine: REAL;
    W: Texts.Writer;

  PROCEDURE Reset*;
  BEGIN q2 := 0.0; q1 := 0.0
  END Reset;

  PROCEDURE Init*;
    VAR k: INTEGER;
      realN, omega: REAL;
  BEGIN realN := FLT(N);
    k := FLOOR(0.5 + ((realN * TargetFrequency) / SamplingRate));
    omega := (2.0 * pi * FLT(k)) / realN;
    sine := Math.sin(omega); cosine := Math.cos(omega);
    coeff := 2.0 * cosine;

    Texts.WriteString(W, "For sampling-rate = "); Texts.WriteReal(W, SamplingRate, 12);
    Texts.WriteString(W, " n = "); Texts.WriteInt(W, N, 4);
    Texts.WriteString(W, " and freq = "); Texts.WriteReal(W, TargetFrequency, 12); Texts.Write(W, ",");
    Texts.WriteLn(W);
    Texts.WriteString(W, "k = "); Texts.WriteInt(W, k, 4); Texts.WriteString(W, " and coeff = ");
    Texts.WriteReal(W, coeff, 12); Texts.WriteLn(W); Texts.WriteLn(W);

    Texts.Append(Oberon.Log, W.buf);
    Reset
  END Init;

  PROCEDURE ProcessSample*(s: Sample);
    VAR q0: REAL;
  BEGIN q0 := coeff * q1 - q2 + FLT(s);
    q2 := q1; q1 := q0
  END ProcessSample;

  PROCEDURE GetRealImag(VAR realPart, imagPart: REAL);
  BEGIN realPart := q1 - q2 * cosine;
    imagPart := q2 * sine
  END GetRealImag;

  PROCEDURE GetMagnitudeSquared*(): REAL;
    VAR result: REAL;
  BEGIN
    result := q1 * q1 + q2 * q2 - q1 * q2 * coeff; RETURN result
  END GetMagnitudeSquared;

BEGIN Texts.OpenWriter(W);
END Goertzel.

MODULE TestGoertzel;
  IMPORT Goertzel, Texts, Math, Oberon;

  VAR W: Texts.Writer;
    testdata: ARRAY Goertzel.N OF Goertzel.Sample;

  PROCEDURE GetStep(f: REAL): REAL;
    VAR result: REAL;
  BEGIN
    result := f * ((2.0 * Goertzel.pi) / Goertzel.SamplingRate); RETURN result
  END GetStep;

  PROCEDURE Generate(f: REAL);
    VAR i: INTEGER;
      step: REAL;
  BEGIN
    step := GetStep(f);
    FOR i := 0 TO Goertzel.N - 1 DO
      testdata[i] := FLOOR(100.0 * Math.sin(FLT(i) * step) + 100.0);
    END
  END Generate;

  PROCEDURE ProcessSamples;
    VAR index: INTEGER;
  BEGIN
    FOR index := 0 TO Goertzel.N - 1 DO
      Goertzel.ProcessSample(testdata[index])
    END
  END ProcessSamples;

  PROCEDURE OptimizedGoertzel;
    VAR magsq, mag: REAL;
  BEGIN
    magsq := Goertzel.GetMagnitudeSquared();
    Texts.WriteString(W, "Relative magnitude squared = "); Texts.WriteReal(W, magsq, 12); Texts.WriteLn(W);
    mag := Math.sqrt(magsq);
    Texts.WriteString(W, "Relative magnitude = "); Texts.WriteReal(W, mag, 12); Texts.WriteLn(W);
  END OptimizedGoertzel;

  PROCEDURE GenerateAndTest(f: REAL);
  BEGIN Texts.WriteString(W, "For test frequency "); Texts.WriteReal(W, f, 12); Texts.Write(W, ":");
    Texts.WriteLn(W);
    Generate(f);
    ProcessSamples;
    OptimizedGoertzel;
    Goertzel.Reset
  END GenerateAndTest;

  PROCEDURE Compute*;
    VAR freq: REAL;
  BEGIN Goertzel.Init;
    GenerateAndTest(Goertzel.TargetFrequency - 250.0);
    GenerateAndTest(Goertzel.TargetFrequency);
    GenerateAndTest(Goertzel.TargetFrequency + 250.0);
    Texts.Append(Oberon.Log, W.buf)
  END Compute;

BEGIN Texts.OpenWriter(W)
END TestGoertzel.