Hi
Ich habe ein Problem bei der Implementierung des Goertzel-Algorithmus
mit Festkommaarithmetik auf einem Cortex M3. Ich messe ein Signal mit
einer Abtastfrequenz von Fs = 8 kHz. Das gesuchte Nutzsignal liegt etwa
bei 2,7 kHz. Die momentane Nutzsignalfrequenz wird über einen Timer
gemessen, so dass der Goertzel für das Messignal abgestimmt wird. Das
Eingangssignal wird nicht gesondert "gefenstert", so dass das
Rechteckfenster hier zum Tragen kommt. Dann gilt ja für die
Frequenzauflösug und damit für die Bandbreite des Goertzel-Filters delta
f = Fs/N und bei N = 2000 -> delta f = 4 Hz. Soweit so gut. Nun ist es
aber so, dass die Ausgangswerte des Systems (die Signalleistung, wie
z.B. hier berechnet:
http://en.literateprograms.org/Goertzel_algorithm_%28C%29) stark
variieren, im bereich von 10^1 bis 10^5. Ich habe einmal die
Referenzfrequenz per Hand im Code absichtlich verschoben und dort waren
die Ausgangswerte größtenteils stabil und klein, so wie man es im
Sperrbereich eines Bandpassfilters erwarten könnte.
Ich weiß nun nicht, ob es zu irgenwelchen Überläufen kommt oder nicht.
Gemessen wird das Signal über ein digitales Mikrofon, das mit 24 Bit
abtastet. Der systemeigene Rauschteppich nimmt aber die untersten 8,5
Bit ein, so dass ich die Eingangssamples um 16 Stellen nach rechts
shiffte, um so zu einer 16.16 Festkommadarstellung zu gelangen. Damit
erhoffe ich mir auch Überläufe bei der Addition/Subtraktion zu
verhindern.
Angehängt ist der komplette Code der aktuellen Fassung. Vielleicht hat
ja jemand eine Idee, wie ich das Problem lösen kann. Wenn ich das Signal
als Rohdaten eine Zeit lang aufnehme und dann in Matlab eine FFT davon
erzeuge, sieht man ganz eindeutig einen starken Peak bei der gesuchten
Frequenz. Der Knackpunkt ist also "nur" der Goertzel Algorithmus auf dem
Mikrocontroller. In Matlab und Floating Point Rechnung liefert der
Goertzel das selbe Ergebnis wie die FFT für die eine Frequenz.
Vielen Dank schon einmal für jede konstruktive Hilfe :-)
Hi Andre,
Die Pole des Goertzel sind auf dem Einheitskreis, deswegen ist das
Problem relativ typisch bei einer geringen Fixkomma-Genauigkeit.
Stichwort "instabil", Energie kommt u.U. einfach mal irgendwoher, je
nach Probe-Frequenz.
Dem entgegensteuern kann man auf verschiedene Weise:
- Pole immer nur auf gewisse quantisierte Stellen setzen
- 'Schlau' runden
Gibt dazu ein paar Papers, Stichworte "Goertzel quantization error"
Ansonsten muss man sich seine Arithmetik scharf anschauen. Bei einer
klassisch in C implementierten Fixkomma-Arithmetik runden die meisten
Implementationen beim Schieben nicht, damit schleppt man also kleine
Quäntchen an energetisch wirksamem Quantisierungsrauschen ein. Bisschen
besser machen das manche DSPs, die spezielle Rundungsoptionen besitzen.
Kann man auch alles in C emulieren, nur wird dann der Algorithmus u.U.
langsam. Bin deswegen dazu übergegangen, die ganzen Filterbänke erst in
DSP Assembler zu prototypen, dann im FPGA zu implementieren, da man da
auch die Arithmetik vollständig unter Kontrolle hat.
Grüsse,
- Strubi
Hi,
was mir auffällt:
- beim cos-Faktor muss Faktor 2 mit rein, machst Du das?
- nullst du v[0] bis v[2] vorher?
- input = (I2SPuffer[((oldIndex) & BUFFER_MASK)][i]>>16);
schmeisst von den 24 Bit 16 weg, nicht gut
- der eine oder andere shift ist sinnfrei (<<2 und danach >>16), die
Umkopiererei kann man auch sparen.
- der Goertzel macht auf jeden Fall nen overflow, das muss dann mit
Deinem N passen.
long long ist bei Euch 64 Bit? Du kannst mit den debug-printf ja auch
komfortabel overflows reporten.
Vllt. einen Skalierungsfaktor mitführen: Du merkst Dir, wie oft Du die
v[0] bis v[2] nach rechts geshiftet hast, um Overflows zu vermeiden.
Cheers
Detlef
Hi und vielen Dank erst einmal für eure Antworten
Martin S. schrieb:> Die Pole des Goertzel sind auf dem Einheitskreis, deswegen ist das> Problem relativ typisch bei einer geringen Fixkomma-Genauigkeit.> Stichwort "instabil", Energie kommt u.U. einfach mal irgendwoher, je> nach Probe-Frequenz.> Dem entgegensteuern kann man auf verschiedene Weise:> - Pole immer nur auf gewisse quantisierte Stellen setzen> - 'Schlau' runden>> Gibt dazu ein paar Papers, Stichworte "Goertzel quantization error"> Ansonsten muss man sich seine Arithmetik scharf anschauen. Bei einer> klassisch in C implementierten Fixkomma-Arithmetik runden die meisten> Implementationen beim Schieben nicht, damit schleppt man also kleine> Quäntchen an energetisch wirksamem Quantisierungsrauschen ein. Bisschen> besser machen das manche DSPs, die spezielle Rundungsoptionen besitzen.> Kann man auch alles in C emulieren, nur wird dann der Algorithmus u.U.> langsam.
Sowas hatte ich mir auch schon überlegt, war allerdings erstaunt, dass
der Goertzel Algorithmus sich aber z.B. mit 8 und 16 Bit Wordbreite
implementieren lässt. TI hatte ja ein Paper zur DTMF Detektion
veröffentlicht, in dem auch nur mit 8 bzw. maximal 16 Bit gearbeitet
wird. Ich bin nun dazu übergegangen im 6.28 Format zu arbeiten. Damit
nehme ich zwar auch den Noise Floor des Systems mit in die Rechnung mit
ein, vermeide dafür aber zusätzliche Probleme durch fehlendes Runden
bzw. Abschneiden der letzten Bits und habe nach oben Kapazitäten gegen
Überläufe.
Der Knackpunkt ist ja, so denke ich es mir, die Wortbreite der
Koeffizienten bzw. des einen, eben 2*cos(2*pi*(f/Fs)). Den belasse ich
im 1.31 Format und muss dann "nur" im Produkt mit dem
Genauigkeitsverlust leben, nicht aber schon bei der Darstellung des
Koeffizienten. Ich erhoffe mir, dass das zu einem besseren Ergebnis
führt. Ich habe aber hier nochmal ein sehr interessantes Paper dazu
gefunden, in dem auf die Quantisierungseffekte beim Goertzel eingegangen
wird:
http://www.adv-radio-sci.net/7/73/2009/ars-7-73-2009.pdfMartin S. schrieb:> Bin deswegen dazu übergegangen, die ganzen Filterbänke erst in> DSP Assembler zu prototypen, dann im FPGA zu implementieren, da man da> auch die Arithmetik vollständig unter Kontrolle hat.
Die Optionen bestehen leider nicht. Die Hardware und alle anderen
Vorgaben sind mittlerweile fix, so dass kein großer Spielraum mehr
existiert.
Detlef _a schrieb:> - beim cos-Faktor muss Faktor 2 mit rein, machst Du das?
Mittlerweile bin ich zu einer Implementierung nach diesem Paper
übergegangen. Daher sind einige Punkte nicht mehr ganz gültig.
http://www.date-conference.com/proceedings/PAPERS/2008/DATE08/PDFFILES/01.5_1.PDFDetlef _a schrieb:> - nullst du v[0] bis v[2] vorher?
Mache ich.
Detlef _a schrieb:> - input = (I2SPuffer[((oldIndex) & BUFFER_MASK)][i]>>16);> schmeisst von den 24 Bit 16 weg, nicht gut
Die Samples werden in 32 Bit Variablen gespeichert. Ein Shift nach
rechts um 16 sollte also die oberen 16 Bits an die untersten
verschieben.
>>> Ein Shift nach rechts um 16 sollte also die oberen 16 Bits an die untersten
verschieben.
Yo, macht er, und die unteren 16Bit sind weg.
Angehängt eine funktionierende Goertzel Routine in fast nur integer,
'winkel' ist der atan2. Nicht schön (eine float-Mult am Ende), nicht
optimal (diese Umkopiererei, die Goertzel-Runde nach der Schleife), aber
gehen tut sie.
Cheers
Detlef
/* ************************************************** */
uint16_t goertzel_21(int16_t r[])
/* ************************************************** */
{
// macht goertzel auf der Frequenz 21/128
int32_t q0,q1,q2;
int16_t k;
q1=0l;
q2=0l;
for(k=0;k<128;k++){
q0=q1+((q1*231)>>13)-q2+r[k];
q2=q1;
q1=q0;
}
q0=q1+((q1*231)>>13)-q2;
q2=q1;
q1=q0;
return (
winkel((int32_t)(q2 *sin(2*M_PI*21/128.0)),
(int32_t)(q1-q2*cos(2*M_PI*21/128.0)))
);
}
Detlef _a schrieb:> Yo, macht er, und die unteren 16Bit sind weg.
Wobei effektiv "nur" 8 Sample Bits dabei verloren gehen. Die untersten 8
Bit sind eh Nullen.
Detlef _a schrieb:> Angehängt eine funktionierende Goertzel Routine in fast nur integer,> 'winkel' ist der atan2. Nicht schön (eine float-Mult am Ende), nicht> optimal (diese Umkopiererei, die Goertzel-Runde nach der Schleife), aber> gehen tut sie.
Ich bräuchte eine Frequenzauflösung von 4 bis 1 Hz, d.h. ich müsste bei
Fs = 8kHz zwischen 2000 und 8000 Abtastwerte verarbeiten. Würde das
deine Implementierung auch hinbekommen? Zumal sie ja eh in in einer 16
Bit Darstellung läuft. Wie stabil ist das Teil denn bei dir? Ich
vermute, dass es bei mir häufig zu überläufen kommt (warum auch immer),
weswegen ich überlege parallel das ganze auf dem STM32F4 Discovery in
float bzw. schlimmstenfalls über Software in double zu implementieren.
Die Frage ist nur, ob sich der Aufwand lohnen würde und ich endlich
stabile Verhältnisse bekomme. Das Problem ist auch, dass sich die
Referenzfrequenz immer etwas ändert (im Bereich von einigen Hz) und ich
somit kein festes Fk angeben kann.
André M. schrieb:> Wie kommt es eigentlich zu q1*231 ?
Bei Dir ist f/fs = 2.7/8. 2*cos(2*pi*f/fs) ist dann -1.044997129, das
ist der Faktor für das q1. -1.044997129 ~ (-1-369/8192) also lautet die
Zeile für Deine Verhältnisse:
q0=-q1-((q1*369)>>13)-q2+r[k];
Overflows: worst case ist der voll ausgesteuerte Sinus auf der Frequenz
des Filters, da muss Du probieren, ob das in 32 Bit reinpasst. Wenn es
knallt, dann bei der Berechnung des Faktors 369/8192, da braucht man 13
zusätzliche Bits wg. 8192. Das kann man aber auch, wie schon gesagt,
monitoren und die q's ggf. eins runtershiften.
>>> Der gemoddete Code ist nicht verifiziert, ich kenne mich und weiß deswegen,
dass da wahrscheinlich noch ein Fehler drin ist ;-(
<<<<
Teste das, bevor Du es benutzt: Einen voll ausgesteuerten Sinus
2.7kHz@8ks/s 8192 samples lang reinfüttern, da muss dasgleiche wie bei
der FFT rauskommen.
Cheers
Detlef
/* ************************************************** */
uint16_t goertzel(int16_t r[])
/* ************************************************** */
/* vllt. buggy ***************************************/
{
// macht goertzel auf der Frequenz 2700/8000
int32_t q0,q1,q2;
int16_t k;
q1=0l;
q2=0l;
for(k=0;k<8192;k++){
q0=-q1-((q1*369)>>13)-q2+r[k];
q2=q1;
q1=q0;
}
q0=-q1-((q1*369)>>13)-q2;
q2=q1;
q1=q0;
return (
winkel((int32_t)(q2 *sin(2*M_PI*2700/8000)),
(int32_t)(q1-q2*cos(2*M_PI*2700/8000)))
);
}
