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.pdf Martin 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.PDF Detlef _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.
> 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; Wie kommt es eigentlich zu q1*231 ?
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))) ); }
Thread beobachten
Seitenaufteilung abschaltenBitte 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
Mit Google-Account einloggen
Noch kein Account? Hier anmelden.