Forum: Mikrocontroller und Digitale Elektronik Goertzel Algorithmus funktioniert nicht


von goertzel (Gast)


Lesenswert?

Hallo

Ich versuche gerade den Goertzel Algorithmus für 440Hz zu 
implementieren. Den Algorithmus habe ich nach dieser Seite 
http://www.answers.com/topic/goertzel-algorithm in C implementiert:
1
while(1)
2
{    
3
    uint16_t coeff = 491; //mit 256 skaliert
4
    int16_t s_prev = 0, s_prev2 = 0;
5
    for(uint16_t i = 0; i < 4096; i++)
6
    {
7
      int16_t sample = get_next_ADC_value();
8
      int16_t s = sample + ((int32_t)(coeff * s_prev) / 256) - s_prev2 ;
9
      s_prev2 = s_prev;
10
      s_prev = s ;
11
    }
12
    s_prev /= 256;
13
    s_prev2 /= 256;
14
    uint16_t power = s_prev2 * s_prev2 + s_prev * s_prev - ((int32_t)coeff * s_prev * s_prev2 / 256);
15
    while (!(UCSR0A & (1<<UDRE0)));
16
    UDR0 = power / 256;
17
}

get_next_ADC_value gibt den nächsten ADC-Wert vorzeichenbehaftet zurück. 
Die ADC läuft im Free-Running Modus und funktioniert einwandfrei, so wie 
der Uart auch.

Allerdings gibt der Algorithmus meistens 0, manchmal auch Werte bis 
leicht über 10 zurück. Allerdings macht es keinen Unterschied ob ich den 
440 Hz Sinus abspiele oder nicht (ADC hängt über einen tiefpass am PC). 
Wenn ich den ADC fest auf 0 lege, ergibt der Algo auch 0. Dass er im 
"normalen" Zustand teilweise leicht über 0 liegt, liegt daran, dass der 
ADC zwischen 0 und -1 schwankt.

von Maxx (Gast)


Lesenswert?

Wie hast du die sample-frequency bestimmt?
(Geht in die 491 der coeff ein)

Die Stützpunkte werden äquidistant vorausgesetzt. Ich erkenne jetzt auf 
dem kurzen Stück (Es fehlen wichtige Teile um das endgültig beurteilen 
zu können) die Dinge, die das sicherstellen.

Ohne diese fehlt dann die Verbindung von Sample-Domäine zur Zeitdomäne 
und damit zu den Frequenzen.

von goertzel (Gast)


Lesenswert?

stimmt, die Samplefrequenz hatte ich vergessen zu erwähnen. Die 
Samplefrequenz liegt bei ~9,6kHz (16Mhz / 128 (ADC-Teiler) / 13 
(Wandlungsdauer)).

Der restliche Code ist nur initialisierung und die get_next_adc_value(). 
Der ADC wird über einen 16-Sample FIFO gepuffert und ist getestet, indem 
ich die Samples über den uart ausgegeben und am PC wieder abgespielt 
habe.

von wutz (Gast)


Lesenswert?

seh ich das richtig, dass du deine Variablen da in einer Schleife 
initialisierst?

sowas macht man doch nur 1 einziges mal?!

von goertzel (Gast)


Lesenswert?

Den coeff könnte man rausnehmen. Die anderen beiden müssen jedes mal auf 
0 gesetzt werden.

von wutz (Gast)


Lesenswert?

goertzel schrieb:
> Den coeff könnte man rausnehmen. Die anderen beiden müssen jedes mal auf
> 0 gesetzt werden.

schon, aber das Initialisieren/Deklarieren einer Variable tut man doch 
nur 1 mal irgendwo am Anfang der Funktion und nicht in einer Schleife.

von nicht Gast (Gast)


Lesenswert?

wird die schleife auch 9600mal in der sekunde durchlaufen???
du musst sicherstellen das get_next_ADC_value() genau alle 104us 
aufgerufen wird....

von Oha (Gast)


Lesenswert?

Ich wuerd den ADC temporaer rausnehmen, und die Daten per UART 
einspielen. So hat man die volle Kontrolle was denn transformiert werden 
soll.

von goertzel (Gast)


Lesenswert?

Damit ihr es mir auch glaubt ;)
1
uint16_t get_next_ADC_value(void)
2
{
3
  while(adc_fifo_read == adc_fifo_write);
4
  uint16_t ret = adc_fifo[adc_fifo_read++];
5
  adc_fifo_read %= ADC_FIFO_LEN;
6
  return 512 - ret;
7
}

Die Samplerate ist konstant, egal wie oft diese Funktion aufgerufen 
wird.

von Maxx (Gast)


Lesenswert?

Ich sehe gerade
1
s_prev /= 256;
2
    s_prev2 /= 256;
3
    uint16_t power = s_prev2 * s_prev2 + s_prev * s_prev - ((int32_t)coeff * s_prev * s_prev2 / 256);

Mit dem ersten beiden /256 verlierst du die Nachkommastellen.
Warum agierst du nicht hier analog in dem du den Datentyp erst 
erweiterst, multiplizierst und die Anzahl der Nachkommstellen danach 
wieder korrigierst, wie du es in den anderen Rechnungen auch machst?

Hier verlierst du im Moment jede Menge Daten.

von goertzel (Gast)


Lesenswert?

Maxx schrieb:
> Hier verlierst du im Moment jede Menge Daten.

Ich denke nicht, dass ich das hier so genau brauche. Ich habe es 
trotzdem getestet - es ändert jedoch nichts am Ergebnis.

von Maxx (Gast)


Lesenswert?

Ich denke schon ;-)

Deine Samples, wenn ich das richtig interpretiere liegen erwartungsgemäß 
zwischen -511 bis 512. Also von -1.996 bis +2.000. Schmeißt du die 
Nachkommstellen weg bleiben nur noch -1,0,1,2 übrig.

Nehmen wir an du rechnest mit der Genauigkeit deiner Samples weiter, 
dann liegen deine Ergebnisse zwischen 0.0 (s_prev, s_prev2 == 0) und 
15,64 (2,-1.996). Du erwartest also bereits keinen besonders hohen 
Werte. Die verbleibenden 4 Bit Auflösung sind die, die du beobachtest.

von goertzel (Gast)


Lesenswert?

Aber diese sollten nicht schwanken wenn nur Werte +-1 am ADC liegen. Und 
sie sollten bei einem 440Hz Sinus nahezu konstant sein, was sie beides 
nicht sind. Der Fehler muss irgendwo anders sein. Ich werde die Division 
trotzdem weglassen.

von Maxx (Gast)


Lesenswert?

Das sind jedenfalls die Restriktionen im Code. Der Goerzel scheint 
ansonsten richtig implementiert.

Offen bleibt dann noch die Frage, ob der adc_fifo auch schnell genug 
ausgelesen und verarbeitet wird. (Was passiert, wenn adc_fifo_write den 
adc_fifo_read einholt?)

Lass ihn doch mal auf einen aufgezeichneten und überprüften Buffer 
laufen (dann hast du Zeit genug um jeden Schritt zu berechnen, ohne dass 
dir der FIFO überläuft/vollläuft)

von Fritz (Gast)


Lesenswert?

Der Goerzel Algorithmus ist vom Prinzip her ein IIR-Filter (siehe deine 
Quellangabe). Laut Wikipedia sind die je nach Implementation sehr 
empfindlich auf Quantisierungsfehler bei Festkommaimplementierung.

goertzel schrieb:
> int16_t sample = get_next_ADC_value();
>       int16_t s = sample + ((int32_t)(coeff * s_prev) / 256) - s_prev2 ;
>       s_prev2 = s_prev;
>       s_prev = s ;

Diese Befehle 4096 mal ausgeführt, kein Wunder das das nur Garbage 
ergibt.

Entweder du implementierst das in float, wenn die CPU schnell genug ist, 
oder du planst die Fixcommaversion besser, wobei ich leicht zweifle ob 
man mit Fixkomma halbwegs genau hinkommt.
Zu Fixkomma: da nimmt man ein Integer und setzt einen (viruellen) 
Dezimalpunkt nach dem höchsten bit (= Vorzeichenbit).
1.0 ergibt sich dann als 0x7fff
So müßtes du deine Daten und Koeffizienten Skalieren. Muß man sich genau 
überlegen in Bezug auf Über- Unterlauf usw.
Das schreit mMn nach Assembler, eine /256 Division ist dann nur ein 
arithmetik-shift-right, üblicherweise nur ein assemblerbefehl, was der 
C-compiler daraus macht ???

von Andreas M. (amesser)


Lesenswert?

Fritz schrieb:

> Das schreit mMn nach Assembler, eine /256 Division ist dann nur ein
> arithmetik-shift-right, üblicherweise nur ein assemblerbefehl, was der
> C-compiler daraus macht ???

Genau das gleiche, auch einen Shift. Compiler sollen schlauer sein als 
man denkt...

Zum Thema: Könnte

uint16_t power = s_prev2 * s_prev2 + s_prev * s_prev - ((int32_t)coeff * 
s_prev * s_prev2 / 256);

einen Überlauf produzieren? Das würde Erklären warum man nur ein 
Rauschen sieht. Vielleicht bin ich doof, kenne den Alogrithmus nicht, 
aber warum wird nochmal mit 256 skaliert? Coeff ist doch auch schon 
skaliert laut Kommentar oben?

generiere Dir doch mal Testdaten für den ADC - d.h. was man mit oder 
ohne Signal erwarten würde und compiliere das Programm für den PC. Da 
kannst Du leichter debuggen und sehen was mit den Variablen passiert.

von goertzel (Gast)


Lesenswert?

Ich habe es jetzt am PC getestet und es läuft einwandfrei. Ich werde 
noch ein bisschen experimentieren wie ich die Rechnung auf Fixkomma 
beschränken kann.

von stm (Gast)


Lesenswert?

> Genau das gleiche, auch einen Shift. Compiler sollen schlauer sein als
> man denkt...

hm bist du dir da sicher? In AVR-Studio mit GCC habe ich mal das Problem 
gehabt, dass er das nicht gemacht hat. Ich musste da die Divisionen 
explizit als Shift-Operationen darstellen!.

Also /256 sollte man durch >>8 ersetzen.

von Andreas M. (amesser)


Lesenswert?

goertzel schrieb
> Ich habe es jetzt am PC getestet und es läuft einwandfrei. Ich werde
> noch ein bisschen experimentieren wie ich die Rechnung auf Fixkomma
> beschränken kann.

Hast du das Programm am PC mit Optimierungen übersetzt? Wenn Du nämlich 
Pech hast, dann macht der PC alle Rechnungen nur in seinen Registern und 
die sind meist 32 Bit breit. Erst am Ende wird auf die geforderte 
Genauigkeit heruntergerechnet. -> Du siehst den Fehler nicht weil er am 
PC nicht auftritt.
Probiere mal alle Variablen am PC als 'volatile' zu kennzeichnen. Dabei 
fällt mir gerade ein, deine Ringpuffer-Variablen sind hoffentlich auch 
alle als 'volatile' markiert?


stm schrieb:
> hm bist du dir da sicher? In AVR-Studio mit GCC habe ich mal das Problem
> gehabt, dass er das nicht gemacht hat. Ich musste da die Divisionen
> explizit als Shift-Operationen darstellen!.

Ja das macht der GCC nicht nur bei AVRs, sondern auch bei ARM, x86...... 
Man muss dazu aber auch die Optimierungen anschalten, '-Os' oder '-O2'. 
Es gibt übrigens noch ganze anderer Optimierungen, die der Compiler 
macht, die würde man als Mensch nur mit viel viel Aufwand hin bekommen - 
z.B. Optimierung der ARM Register-Interlocks.

von goertzel (Gast)


Lesenswert?

Ich habe die das Programm erst mit double geschrieben und das dann durch 
interger verschiedener Größen ersetzt. Bei 16bit funktioniert es bei 
wenigen Samples noch gut, da ich aber ein kleineres Fenster brauche, 
werde ich auf 24bit/32bit gehen.

PS: Der Unterschied zwischen 16 und 32bit war bemerkbar, da wurde nichts 
wegoptimiert.

von goertzel (Gast)


Lesenswert?

Es funktioniert jetzt auch mit mehreren Frequenzen gleichzeitig, 
allerdings schleichen sich ab mehr als 2 Frequenzen Fehler ein, da der 
FIFO Buffer nicht mehr reicht, bzw. der C-Code zu lahm ist. Ich werde 
das ganze wohl doch in Assembler implementieren.

von Johann L. (gjlayde) Benutzerseite


Lesenswert?

Andreas Messer schrieb:
> Fritz schrieb:
>
>> Das schreit mMn nach Assembler, eine /256 Division ist dann nur ein
>> arithmetik-shift-right, üblicherweise nur ein assemblerbefehl, was der
>> C-compiler daraus macht ???
>
> Genau das gleiche, auch einen Shift. Compiler sollen schlauer sein als
> man denkt...

Tatsächlich ist der Compiler schlauer, als ihr denkt. Denn:

Die Division von oben isk kein Shift um 8 nach rechts:
1
int div256 (int a)
2
{
3
    return a / 256;
4
}
 
wird mit avr-gcc zu
 
1
div256:
2
  sbrs r25,7
3
  rjmp .L2
4
  subi r24,1
5
  sbci r25,-1
6
.L2:
7
  mov r24,r25
8
  lsl r25
9
  sbc r25,r25
10
  ret
 
Der Code ist korrekt und geht auch nicht kürzer; zumindest seht ich auf 
die Schnelle nicht, wo da noch was gespart werden kann.

Und wenn... immer noch besser suboptimaler Code als falscher Code ;-)

Bitte melde dich an um einen Beitrag zu schreiben. Anmeldung ist kostenlos und dauert nur eine Minute.
Bestehender Account
Schon ein Account bei Google/GoogleMail? Keine Anmeldung erforderlich!
Mit Google-Account einloggen
Noch kein Account? Hier anmelden.