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:
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.
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.
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.
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.
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.
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.
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.
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.
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)
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 ???
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.
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.
> 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.
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.
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.
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.
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
intdiv256(inta)
2
{
3
returna/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 ;-)