Goertzel Algorithmus

[Til H. (turakar)]

Hallo,
ich habe den Goertzel Algorithmus in Java implementiert, um ihn zu 
testen. Nun ist das Problem, dass ich nicht 100% sicher bin, wie ich die 
Magnituden zu verstehen habe. Es geht grundsätzlich darum, zu erkennen, 
ob ein Ton in dem Signal vorkommt oder nicht (eine gewisse Hemmschwelle 
muss natürlich vorhanden sein). Auch bin ich mir nicht ganz sicher, wie 
ich den Debug-Ton zu generieren habe, jedoch wirkt der im Debugger schon 
ganz richtig. Bei meiner Implementierung generiere ich nun künstlich 
einen Ton, um ihn dann mit dem Goertzel auf verschiedenen Frequenzen 
abzutasten. Nun erwarte ich, dass ich die höchste Magnitude bekomme, 
wenn ich nach der gleichen Frequenz suche (941Hz). Allerdings liegt die 
höchste Amplitude bei 961Hz, was mir nach ein bisschen viel Abweichung 
aussieht, oder ist das im akzeptablen Rahmen? (Ein double hat in Java 
64-bit Genauigkeit)

package de.meinen_nachnamen_bekommt_ihr_nicht.test.goertzel;

public class TestGoertzel {
  
  /**
   * constant for 2 * pi
   */
  private static final double TWO_PI = 2d * Math.PI;
  
  /**
   * the k factor, precomputed
   */
  private int k;
  /**
   * the sine of omega, precomputed
   */
  private double sin;
  /**
   * the cosine of omega, precomputed
   */
  private double cos;
  /**
   * cos * 2, precomputed
   */
  private double coeff;
  
  /**
   * temporarily storage variables
   */
  private double q0, q1, q2;
  
  /**
   * default constructor, does nothing
   */
  public TestGoertzel() {
    
  }
  
  /**
   * sets the constants to work with NOTE: targetFrequency % (blocksize / samplingRate) should be 0
   * @param targetFrequency the frequency to listen to
   * @param blocksize the amount of measurements in one cycle
   * @param samplingRate the sampling rate of your measurements
   */
  public void setConstants(int targetFrequency, int blocksize, int samplingRate) {
    k = (int) (0.5f + (blocksize * targetFrequency) / ((double) samplingRate));
    double omega = (TWO_PI / blocksize) * k;
    sin = Math.sin(omega);
    cos = Math.cos(omega);
    coeff = 2 * cos;
  }
  
  /**
   * starts a measurement
   */
  public void start() {
    q1 = 0;
    q2 = 0;
  }
  
  /**
   * processes one sample
   * @param sample the sample
   */
  public void processSample(double sample) {
    q0 = coeff * q0 - q2 + sample;
    q2 = q1;
    q1 = q0;
  }
  
  /**
   * ends a measurement and computes real and imaginary magnitude
   * @param if null or smaller than 2, a new array is used
   * @return result [0] is real and [1] is imaginary
   */
  public double[] end(double[] results) {
    if(results == null || results.length < 2) {
      results = new double[2];
    }
    results[0] = q1 - q2 * cos;
    results[1] = q2 * sin;
    return results;
  }
  
  /**
   * computes the magnitude based on real and imaginary magnitude
   * @param realImag the real and imag array given by end()
   * @return magnitude
   */
  public double computeMagnitude(double[] realImag) {
    return computeMagnitude(computeMagnitudeSquared(realImag));
  }
  
  /**
   * computes the magnitude based on the squared magnitude (Same as Math.sqrt(magSquared))
   * @param magSquared the magnitude^2
   * @return magnitude
   */
  public double computeMagnitude(double magSquared) {
    return Math.sqrt(magSquared);
  }
  
  /**
   * computes the squared magnitude based on the real and imaginary magnitudes
   * @param realImag
   * @return magnitude^2
   */
  public double computeMagnitudeSquared(double[] realImag) {
    return realImag[0] * realImag[0] + realImag[1] * realImag[1];
  }

  /**
   * for testing purposes
   * @param args
   */
  public static void main(String[] args) {
    int blocksize = 205;
    int freq = 941;
    double amplitude = 1;
    int samplingRate = 8000;
    
    double[] samples = new double[blocksize];
    for(int i = 0; i < blocksize; i++) {
        samples[i] = Math.sin(2 * Math.PI * freq * i * (1f / samplingRate)) * amplitude;
    }
    
    System.out.println("Frequency: " + freq + " Amplitude: " + amplitude + " Sampling Rate: " + samplingRate + " Blocksize: " + blocksize);
    int testFreqMin = 641, testFreqMax = 1241, testFreqStep = 5;
    double highestMagnitude = 0, highestMagnitudeFreq = 0;
    TestGoertzel goertzel = new TestGoertzel();
    for(int testFreq = testFreqMin; testFreq <= testFreqMax; testFreq += testFreqStep) {
      goertzel.setConstants(testFreq, blocksize, samplingRate);
      goertzel.start();
      for(int i = 0; i < blocksize; i++) {
          goertzel.processSample(samples[i]);
      }
      double[] results = goertzel.end(null);
      double mag = goertzel.computeMagnitude(results);
      if(highestMagnitude < mag) {
        highestMagnitudeFreq = testFreq;
        highestMagnitude = mag;
      }
      System.out.println("testFreq " + testFreq + ": " + mag);
    }
    System.out.println("Highest Magnitude by frequency: " + highestMagnitudeFreq);
  }

}

Output:

Frequency: 941 Amplitude: 1.0 Sampling Rate: 8000 Blocksize: 205
testFreq 641: 1.1014216897838298
testFreq 646: 3.8596151492613955
testFreq 651: 4.29501571613194
testFreq 656: 4.295015716131877
testFreq 661: 4.29501571613189
testFreq 666: 4.295015716131897
testFreq 671: 4.295015716131916
testFreq 676: 4.295015716131925
testFreq 681: 4.295015716131909
testFreq 686: 4.218815832793116
testFreq 691: 4.786483369227722
testFreq 696: 4.7864833692276685
testFreq 701: 4.786483369227665
testFreq 706: 4.786483369227664
testFreq 711: 4.786483369227692
testFreq 716: 4.786483369227652
testFreq 721: 4.786483369227655
testFreq 726: 4.705740156350941
testFreq 731: 5.477180719536364
testFreq 736: 5.4771807195363165
testFreq 741: 5.477180719536358
testFreq 746: 5.477180719536329
testFreq 751: 5.477180719536303
testFreq 756: 5.477180719536324
testFreq 761: 5.399417600903364
testFreq 766: 6.512243587134702
testFreq 771: 6.512243587134612
testFreq 776: 6.5122435871346385
testFreq 781: 6.5122435871346385
testFreq 786: 6.5122435871346385
testFreq 791: 6.5122435871346385
testFreq 796: 6.5122435871346385
testFreq 801: 6.464394129158938
testFreq 806: 8.22339474779077
testFreq 811: 8.22339474779065
testFreq 816: 8.223394747790651
testFreq 821: 8.22339474779065
testFreq 826: 8.223394747790667
testFreq 831: 8.223394747790628
testFreq 836: 8.223394747790643
testFreq 841: 8.322212220019779
testFreq 846: 11.57044305758967
testFreq 851: 11.57044305758951
testFreq 856: 11.570443057589522
testFreq 861: 11.570443057589484
testFreq 866: 11.570443057589502
testFreq 871: 11.570443057589532
testFreq 876: 11.570443057589522
testFreq 881: 12.673570150608604
testFreq 886: 20.962568122599787
testFreq 891: 20.962568122599706
testFreq 896: 20.962568122599592
testFreq 901: 20.962568122599635
testFreq 906: 20.96256812259962
testFreq 911: 20.96256812259965
testFreq 916: 20.96256812259962
testFreq 921: 94.82873021927618
testFreq 926: 196.71534868622342
testFreq 931: 196.71534868621714
testFreq 936: 196.7153486862186
testFreq 941: 196.71534868621825
testFreq 946: 196.71534868621708
testFreq 951: 196.71534868621723
testFreq 956: 196.71534868621723
testFreq 961: 206.73680252933343
testFreq 966: 23.91530522974619
testFreq 971: 23.915305229741293
testFreq 976: 23.915305229741303
testFreq 981: 23.91530522974124
testFreq 986: 23.915305229741303
testFreq 991: 23.91530522974124
testFreq 996: 22.93319093242113
testFreq 1001: 10.707712132182861
testFreq 1006: 10.707712132183282
testFreq 1011: 10.70771213218326
testFreq 1016: 10.707712132183273
testFreq 1021: 10.70771213218326
testFreq 1026: 10.707712132183268
testFreq 1031: 10.707712132183271
testFreq 1036: 10.118297166617303
testFreq 1041: 6.664075687565079
testFreq 1046: 6.664075687565156
testFreq 1051: 6.664075687565142
testFreq 1056: 6.664075687565144
testFreq 1061: 6.664075687565142
testFreq 1066: 6.664075687565144
testFreq 1071: 6.664075687565142
testFreq 1076: 6.26097813381634
testFreq 1081: 4.7120335261679935
testFreq 1086: 4.712033526168036
testFreq 1091: 4.7120335261680495
testFreq 1096: 4.712033526168047
testFreq 1101: 4.712033526168042
testFreq 1106: 4.7120335261680495
testFreq 1111: 4.712033526168047
testFreq 1116: 4.409963532393034
testFreq 1121: 3.568203590880409
testFreq 1126: 3.5682035908804295
testFreq 1131: 3.568203590880427
testFreq 1136: 3.568203590880427
testFreq 1141: 3.568203590880427
testFreq 1146: 3.568203590880427
testFreq 1151: 3.568203590880427
testFreq 1156: 3.3284704579910462
testFreq 1161: 2.821336328207954
testFreq 1166: 2.8213363282079578
testFreq 1171: 2.8213363282079547
testFreq 1176: 2.8213363282079533
testFreq 1181: 2.8213363282079547
testFreq 1186: 2.8213363282079533
testFreq 1191: 2.6238123039485988
testFreq 1196: 2.2991334444623703
testFreq 1201: 2.2991334444623734
testFreq 1206: 2.2991334444623783
testFreq 1211: 2.299133444462382
testFreq 1216: 2.29913344446238
testFreq 1221: 2.299133444462381
testFreq 1226: 2.2991334444623734
testFreq 1231: 2.132255407129717
testFreq 1236: 1.9167604494689559
testFreq 1241: 1.9167604494689596
Highest Magnitude by frequency: 961.0

MfG,
Tilman H.

PS: Gibt es hier Spoiler? Wären hilfreich.

EDIT: Nach ein paar Versuchsreihen stellt sich heraus, dass die 
Abweichung wohl immer bei 20Hz liegt. Ich würde die jetzt jedoch ungerne 
einfach immer abziehen, sondern lieber wissen, wo der Fehler liegt.

[Michael H. (dowjones)]

Sieht aus als fehlte in start() ein q0=0.

[Mr. Tom (Gast)]

Til Hoff schrieb:
> ..., sondern lieber wissen, wo der Fehler liegt.

Hier seitenweise Quellcode und Output zu posten. Wofür gibt es Anhänge. 
Für den Output wäre statt der fast endlosen Liste ein Graph der 
Magnitude über die Frequenz x-mal übersichtlicher.

Dann wäre dir auch aufgefallen, das deine Maximumlage bei 961 Hz nur ein 
Artefakt der Rechnung ist.

[m.u. (Gast)]

Til Hoff schrieb:
> PS: Gibt es hier Spoiler?
Was für Teile?
http://de.wikipedia.org/wiki/Spoiler

[Werner M. (Gast)]

Til Hoff schrieb:
> Nach ein paar Versuchsreihen stellt sich heraus, dass die
> Abweichung wohl immer bei 20Hz liegt.

Lass vor der Suche nach dem Maximum mal einen phasenneutralen 
FIR-Tiefpass (z.B. Gleitenden-Mittelwert mit vielleicht 9 Stützstellen) 
über den Output laufen. Dann sieht die Welt wieder ganz anders aus.
Wen interessieren eigentlich die zwölfunddrölfzigste Nachkommastellen? 
Hast du dazu mal eine Fehlerbetrachtung angestellt und bestimmt, 
wieviele davon relevant sind?

[Til H. (turakar)]

Ich stufe das Problem als eine ungenaue Berechnung ein. (Liebes, liebes, 
Pi)
Im Moment habe ich jedoch eine Frage zum Verhalten des Algorithmuses 
selber, und zwar kommt es bei meiner Implementierung bei US Frequenzen 
zu so hohen Werten in der Berechnung, dass es zu einem Überlauf kommt. 
Wie groß werden denn solche Werte bei 20kHz oä.?

Edit: Beim nächsten mal verwende ich Anhänge.

[Til H. (turakar)]

EDIT: Mein Fehler. Habe auf 200kHz getestet, eine Null zu viel.