Goertzel Algorithmus

Til H. #3479895
Lesenswert?
• ▲
▼ 
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)
1
package de.meinen_nachnamen_bekommt_ihr_nicht.test.goertzel;

2


3
public class TestGoertzel {

4
  

5
  /**

6
   * constant for 2 * pi

7
   */

8
  private static final double TWO_PI = 2d * Math.PI;

9
  

10
  /**

11
   * the k factor, precomputed

12
   */

13
  private int k;

14
  /**

15
   * the sine of omega, precomputed

16
   */

17
  private double sin;

18
  /**

19
   * the cosine of omega, precomputed

20
   */

21
  private double cos;

22
  /**

23
   * cos * 2, precomputed

24
   */

25
  private double coeff;

26
  

27
  /**

28
   * temporarily storage variables

29
   */

30
  private double q0, q1, q2;

31
  

32
  /**

33
   * default constructor, does nothing

34
   */

35
  public TestGoertzel() {

36
    

37
  }

38
  

39
  /**

40
   * sets the constants to work with NOTE: targetFrequency % (blocksize / samplingRate) should be 0

41
   * @param targetFrequency the frequency to listen to

42
   * @param blocksize the amount of measurements in one cycle

43
   * @param samplingRate the sampling rate of your measurements

44
   */

45
  public void setConstants(int targetFrequency, int blocksize, int samplingRate) {

46
    k = (int) (0.5f + (blocksize * targetFrequency) / ((double) samplingRate));

47
    double omega = (TWO_PI / blocksize) * k;

48
    sin = Math.sin(omega);

49
    cos = Math.cos(omega);

50
    coeff = 2 * cos;

51
  }

52
  

53
  /**

54
   * starts a measurement

55
   */

56
  public void start() {

57
    q1 = 0;

58
    q2 = 0;

59
  }

60
  

61
  /**

62
   * processes one sample

63
   * @param sample the sample

64
   */

65
  public void processSample(double sample) {

66
    q0 = coeff * q0 - q2 + sample;

67
    q2 = q1;

68
    q1 = q0;

69
  }

70
  

71
  /**

72
   * ends a measurement and computes real and imaginary magnitude

73
   * @param if null or smaller than 2, a new array is used

74
   * @return result [0] is real and [1] is imaginary

75
   */

76
  public double[] end(double[] results) {

77
    if(results == null || results.length < 2) {

78
      results = new double[2];

79
    }

80
    results[0] = q1 - q2 * cos;

81
    results[1] = q2 * sin;

82
    return results;

83
  }

84
  

85
  /**

86
   * computes the magnitude based on real and imaginary magnitude

87
   * @param realImag the real and imag array given by end()

88
   * @return magnitude

89
   */

90
  public double computeMagnitude(double[] realImag) {

91
    return computeMagnitude(computeMagnitudeSquared(realImag));

92
  }

93
  

94
  /**

95
   * computes the magnitude based on the squared magnitude (Same as Math.sqrt(magSquared))

96
   * @param magSquared the magnitude^2

97
   * @return magnitude

98
   */

99
  public double computeMagnitude(double magSquared) {

100
    return Math.sqrt(magSquared);

101
  }

102
  

103
  /**

104
   * computes the squared magnitude based on the real and imaginary magnitudes

105
   * @param realImag

106
   * @return magnitude^2

107
   */

108
  public double computeMagnitudeSquared(double[] realImag) {

109
    return realImag[0] * realImag[0] + realImag[1] * realImag[1];

110
  }

111


112
  /**

113
   * for testing purposes

114
   * @param args

115
   */

116
  public static void main(String[] args) {

117
    int blocksize = 205;

118
    int freq = 941;

119
    double amplitude = 1;

120
    int samplingRate = 8000;

121
    

122
    double[] samples = new double[blocksize];

123
    for(int i = 0; i < blocksize; i++) {

124
        samples[i] = Math.sin(2 * Math.PI * freq * i * (1f / samplingRate)) * amplitude;

125
    }

126
    

127
    System.out.println("Frequency: " + freq + " Amplitude: " + amplitude + " Sampling Rate: " + samplingRate + " Blocksize: " + blocksize);

128
    int testFreqMin = 641, testFreqMax = 1241, testFreqStep = 5;

129
    double highestMagnitude = 0, highestMagnitudeFreq = 0;

130
    TestGoertzel goertzel = new TestGoertzel();

131
    for(int testFreq = testFreqMin; testFreq <= testFreqMax; testFreq += testFreqStep) {

132
      goertzel.setConstants(testFreq, blocksize, samplingRate);

133
      goertzel.start();

134
      for(int i = 0; i < blocksize; i++) {

135
          goertzel.processSample(samples[i]);

136
      }

137
      double[] results = goertzel.end(null);

138
      double mag = goertzel.computeMagnitude(results);

139
      if(highestMagnitude < mag) {

140
        highestMagnitudeFreq = testFreq;

141
        highestMagnitude = mag;

142
      }

143
      System.out.println("testFreq " + testFreq + ": " + mag);

144
    }

145
    System.out.println("Highest Magnitude by frequency: " + highestMagnitudeFreq);

146
  }

147


148
}

Output:
1
Frequency: 941 Amplitude: 1.0 Sampling Rate: 8000 Blocksize: 205

2
testFreq 641: 1.1014216897838298

3
testFreq 646: 3.8596151492613955

4
testFreq 651: 4.29501571613194

5
testFreq 656: 4.295015716131877

6
testFreq 661: 4.29501571613189

7
testFreq 666: 4.295015716131897

8
testFreq 671: 4.295015716131916

9
testFreq 676: 4.295015716131925

10
testFreq 681: 4.295015716131909

11
testFreq 686: 4.218815832793116

12
testFreq 691: 4.786483369227722

13
testFreq 696: 4.7864833692276685

14
testFreq 701: 4.786483369227665

15
testFreq 706: 4.786483369227664

16
testFreq 711: 4.786483369227692

17
testFreq 716: 4.786483369227652

18
testFreq 721: 4.786483369227655

19
testFreq 726: 4.705740156350941

20
testFreq 731: 5.477180719536364

21
testFreq 736: 5.4771807195363165

22
testFreq 741: 5.477180719536358

23
testFreq 746: 5.477180719536329

24
testFreq 751: 5.477180719536303

25
testFreq 756: 5.477180719536324

26
testFreq 761: 5.399417600903364

27
testFreq 766: 6.512243587134702

28
testFreq 771: 6.512243587134612

29
testFreq 776: 6.5122435871346385

30
testFreq 781: 6.5122435871346385

31
testFreq 786: 6.5122435871346385

32
testFreq 791: 6.5122435871346385

33
testFreq 796: 6.5122435871346385

34
testFreq 801: 6.464394129158938

35
testFreq 806: 8.22339474779077

36
testFreq 811: 8.22339474779065

37
testFreq 816: 8.223394747790651

38
testFreq 821: 8.22339474779065

39
testFreq 826: 8.223394747790667

40
testFreq 831: 8.223394747790628

41
testFreq 836: 8.223394747790643

42
testFreq 841: 8.322212220019779

43
testFreq 846: 11.57044305758967

44
testFreq 851: 11.57044305758951

45
testFreq 856: 11.570443057589522

46
testFreq 861: 11.570443057589484

47
testFreq 866: 11.570443057589502

48
testFreq 871: 11.570443057589532

49
testFreq 876: 11.570443057589522

50
testFreq 881: 12.673570150608604

51
testFreq 886: 20.962568122599787

52
testFreq 891: 20.962568122599706

53
testFreq 896: 20.962568122599592

54
testFreq 901: 20.962568122599635

55
testFreq 906: 20.96256812259962

56
testFreq 911: 20.96256812259965

57
testFreq 916: 20.96256812259962

58
testFreq 921: 94.82873021927618

59
testFreq 926: 196.71534868622342

60
testFreq 931: 196.71534868621714

61
testFreq 936: 196.7153486862186

62
testFreq 941: 196.71534868621825

63
testFreq 946: 196.71534868621708

64
testFreq 951: 196.71534868621723

65
testFreq 956: 196.71534868621723

66
testFreq 961: 206.73680252933343

67
testFreq 966: 23.91530522974619

68
testFreq 971: 23.915305229741293

69
testFreq 976: 23.915305229741303

70
testFreq 981: 23.91530522974124

71
testFreq 986: 23.915305229741303

72
testFreq 991: 23.91530522974124

73
testFreq 996: 22.93319093242113

74
testFreq 1001: 10.707712132182861

75
testFreq 1006: 10.707712132183282

76
testFreq 1011: 10.70771213218326

77
testFreq 1016: 10.707712132183273

78
testFreq 1021: 10.70771213218326

79
testFreq 1026: 10.707712132183268

80
testFreq 1031: 10.707712132183271

81
testFreq 1036: 10.118297166617303

82
testFreq 1041: 6.664075687565079

83
testFreq 1046: 6.664075687565156

84
testFreq 1051: 6.664075687565142

85
testFreq 1056: 6.664075687565144

86
testFreq 1061: 6.664075687565142

87
testFreq 1066: 6.664075687565144

88
testFreq 1071: 6.664075687565142

89
testFreq 1076: 6.26097813381634

90
testFreq 1081: 4.7120335261679935

91
testFreq 1086: 4.712033526168036

92
testFreq 1091: 4.7120335261680495

93
testFreq 1096: 4.712033526168047

94
testFreq 1101: 4.712033526168042

95
testFreq 1106: 4.7120335261680495

96
testFreq 1111: 4.712033526168047

97
testFreq 1116: 4.409963532393034

98
testFreq 1121: 3.568203590880409

99
testFreq 1126: 3.5682035908804295

100
testFreq 1131: 3.568203590880427

101
testFreq 1136: 3.568203590880427

102
testFreq 1141: 3.568203590880427

103
testFreq 1146: 3.568203590880427

104
testFreq 1151: 3.568203590880427

105
testFreq 1156: 3.3284704579910462

106
testFreq 1161: 2.821336328207954

107
testFreq 1166: 2.8213363282079578

108
testFreq 1171: 2.8213363282079547

109
testFreq 1176: 2.8213363282079533

110
testFreq 1181: 2.8213363282079547

111
testFreq 1186: 2.8213363282079533

112
testFreq 1191: 2.6238123039485988

113
testFreq 1196: 2.2991334444623703

114
testFreq 1201: 2.2991334444623734

115
testFreq 1206: 2.2991334444623783

116
testFreq 1211: 2.299133444462382

117
testFreq 1216: 2.29913344446238

118
testFreq 1221: 2.299133444462381

119
testFreq 1226: 2.2991334444623734

120
testFreq 1231: 2.132255407129717

121
testFreq 1236: 1.9167604494689559

122
testFreq 1241: 1.9167604494689596

123
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.
Antwort Zitieren Bearbeiten Thread verschieben Thread sperren Anmeldepflicht aktivieren Anpinnen Thread löschen Thread mit anderem zusammenführen Beitrag melden
Mr. Tom Gast #3482535
Lesenswert?
• ▲
▼ 
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.
Angehängte Dateien:
Antwort Zitieren Bearbeiten Löschen Beitrag melden
Werner M. Gast #3483294
Lesenswert?
• ▲
▼ 
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?
Antwort Zitieren Bearbeiten Löschen Beitrag melden
Til H. #3485165
Lesenswert?
• ▲
▼ 
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.
Antwort Zitieren Bearbeiten Löschen Beitrag melden
Thread beobachten

Antwort schreiben

Bitte melde dich an, um einen Beitrag zu schreiben.

Bestehender Account

oder

Mit Google-Account einloggen

Noch kein Account?

Die Registrierung ist kostenlos und dauert nur eine Minute.

Jetzt registrieren