Goertzel Algorithmus (Frequenz + Phase)

[stoxx]

Hallo liebe Delphianer :)


wisst ihr zufällig, wo man im Goertzel Algorithmus die Phase (und Amplitude) herbekommt?

hab den Algorithmus mal in Delphi umgesetzt, die Frequenzbestimmung funktioniert soweit wie gewünscht. (Projekt im Anhang)
In diesem PDF steht leider nur, dass man die Phase bei der reinen Frequenzbestimmung ja nicht bräuchte .. hmm :(

http://www.ti.com/lit/an/spra096a/spra096a.pdf

ich bin jetzt in Mathe so tief leider nicht gleich so fit, das zu verstehen.
Vielleicht kann ja jemand helfen ?

Delphi-Quellcode:

function Goertzel(const x : TDoubleArray; frequency : Double; SampleRate : Integer) : Extended;
var Skn, Skn1, Skn2 : Extended;
    i : Integer;
    WNk : Extended;
    coeff : Extended;
    power : Extended;
begin

  Skn := 0; Skn1 := 0; skn2 := 0;
  coeff := 2 * Cos( 2 * Pi * frequency / SampleRate);
  for I := 0 to Length(x)- 1  do begin

    Skn :=  x[i] + (coeff * Skn1) - Skn2 ;
    Skn2 := Skn1;
    Skn1 := Skn;

  end; // for I 
 result := Skn2 * Skn2 + Skn1 * Skn1 - coeff * Skn1*skn2;
end;

Danke ! :)

[Peter1999]

Hallo,

guck mal hier,

da steht neben der "schnellen" Berechnung der Amplitude auch die Berechnung der komplexen Transformierten.
Damit hast dann Real- und Imaginärteil, um deine Phase zu berechnen (arctan(Im/Re)).

Viele Grüße...

[stoxx]

Das klappt leider alles immer noch nicht so richtig.
Die Winkelberechnung funktioniert manchmal, wenn die Soundlänge genau 1000 ms lang ist.
hast Du den Algorithmus schonmal umgesetzt?
weißt Du, wo und wie da korrekt normiert werden muss?

[stoxx]

....

[Peter1999]

So,
du hattest mit der Normierung recht, das ist auch leicht unübersichtlich, wenn man sich so im Internet umsieht. Da ich mich sonst auch eher mit FFTs beschäftige, war das jetzt auch etwas ungewohnt und ich übernehme keine Garantie ;) . Egal...

Dass die Phase bei einem einzelnen Signal keine Aussage besitzt, ist dir aber schon klar? Die wird erst sinnvoll, wenn du Frequenzgänge bestimmen willst (mind. 2 synchrone Signale).

Ich habe mal ein eigenes kleines Beispielprogramm geschrieben, bei dem ich nun davon ausgehe, dass es (ungefähr) funktioniert.
Das ganze sollte zwar noch etwas aufgeräumt und ordentlich in eine Klasse verpackt werden, macht aber ansonsten, was es soll.

Auch die Frequenzfindung kann noch deutlich verbessert werden, wenn man weiß, dass man wirklich nur eine diskrete Spektrallinie im Signal hat. Bei einem kontinuierlichen Spektrum kannst du das natürlich vergessen. Dann gibt es nicht viel zu Optimieren (außer der Frequenzauflösung vielleicht?).

Du kannst es dir ja mal anschauen.

Viele Grüße...

[stoxx]

vielen Dank ! .. also erstmal .. vor allen Dingen für die Mühe, ist recht praktisch gemacht, das Programm :-)
die Signalfindung soll nich verbessert werden, es geht um das prüfen einiger (an der hand abzählbarer) Frequenzen.

wie ich jetzt so gelesen habe, soll gerade der Vorteil vom Goertzel Algorithmus sein, dass man an kein festes Fenster gebunden ist, sondern sogar ein beliebiges großes (einziges) Datenarray berechnen kann. (Somit unterscheidet sich das Fenster hier glaub ich etwas)

aber nun zum den Problemen .. siehe bild ..
vor dem Problem steh ich schon die ganze Zeit .. man variiert mal bissl die Samplefrequenz, und die Abtastfrequenz und dann kommt nur Quatsch raus.
Die Amplitude ist ständig falsch (manchmal sogar in Anhängigkeit von der Phase kommt was anderes raus.. )
Die Phase sollte aber auch bei Null beginnen, wie meinst Du das, sie hätte keine Bedeutung? Lässt sie sich nicht korrekt berechnen?

Es scheint, als wäre das alles komplizierter, als gedacht ...
hier wird von Korrekturfaktoren geredet, aber ich versteh nix !

http://www.dsprelated.com/showmessage/87369/1.php

http://www.nntpnews.info/threads/695...ctral-response



Von Mathworks hab ich mal noch in Delphi implementiert, funktioniert aber auch nicht wie gewünscht. Manchmal kommt auch murx raus.
Zum Finden der Frequenzen sind die Verfahren alle geeignent, nur beim richtigen Imaginär und Realteil versagen sie. ich vermute ja immer noch, dass etwas eine Korrektur notwendig ist, oder noch was normiert werden muss, er also irgendwie zu hoch oder zu tief rechnet.

(vielleicht hab ich auch noch einen Fehler drin ... hmm)
ist im anhang, da ist auch ein Algo drin, wenn man die fenster overlappen möchte.
Das wird an anderer Stelle irgendwo empfohlen, da an den Fensterrändern Fehler entstehen ..


http://www.mathworks.com/matlabcentr...ent/goertzel.m

[stoxx]

ahso .. das wichtigste vergessen.

http://www.mstarlabs.com/dsp/goertzel/goertzel.html

dieser Algo im Link unterscheidet sich von dem in Deinem Link etwas.
Der liefert aber wenigistens stabile (also von der reingesteckten Phase) unabhängige Amplituden. (sollte auch so sein :)

da die Formel leicht anders aussieht, vor allen Dingen beim Bestimmen des Im und Re Teils, könnte also auch viel Murks im Internet unterwegs sein...

Allerdings stimmt hier die Phase auch noch nicht !

[Peter1999]

Das Ergebnis im Bild wird doch (höchst wahrscheinlich sogar) stimmen.
Die ermittelte Frequenz liegt bei 193,6 Hz, die tatsächliche aber bei 200 Hz. Damit siehst du nicht die exakte Spektrallinie, sondern ein "verschmiertes" Nebenband. Dass dort die Amplitude deutlich geringer ist, ist ja bei einem einfachen Sinussignal normal. Wo soll die Signalleistung abseits der 200Hz auch herkommen. Nur wenn du die Frequenz besser bestimmst (also möglichst wirklich die 200Hz), bekommst du das exakte Ergebnis.

Gib mal für die Goertzel-Frequenz von Hand die 200Hz vor und klicke auf den Goertzel-Button, dann siehst du, dass es wirklich nur daran liegt, dass du auf dem Seitenband gelandet bist.

Zur besseren Bestimmung der Frequenz musst du aber ein Optimierungsproblem lösen. Das ist es, was ich mit verbesserter Frequenzsuche meinte.

Einfachstes Verfahren:
Der Algorithmus müsste also die 2 größten Amplituden finden. Den Bereich dazwischen wieder in äquidistante (im Frequenzbereich) Stücke einteilen --> da wieder die 2 größten Amplituden finden --> dann wieder neu einteilen usw. ...

Das ganze musst du dann machen bis ein Konvergenzkriterium erreicht ist, z.B. sich die Max.-Amplitude nicht mehr als 0,1% von der letzten Max.-Amplitude unterscheidet.

Wie gesagt, bei der Amplitude ist das durchaus sinnvoll. Die Phase ist aber immer nur die Phase eines Sinus bezüglich deines gewählten Zeitfensters. Was du also damit machen willst, ist mir noch schleierhaft, solange du nicht mehrere Kanäle oder synchrone Daten vergleichen willst.

Du schreibst von einer Korrektur, die nötig sein könnte. Das ist falsch! Die DFT, Goertzel oder FFT sind immer eindeutig und exakt umkehrbar! Was mit der Korrektur gemeint ist, ist folgendes: Man kann aus der Breite des verschmierten Spektrums den Maximalwert abschätzen (siehe Modalanalyse). Man kann auch durch geeignete Fensterfunktionen (Hanning, Blackman,...) das Verschmieren reduzieren. Das würde ich dir aber alles nicht empfehlen, sondern (wie gesagt) lieber die Frequenz richtiger bestimmen. Dann gibt es das Problem nämlich nicht.

Zusammenfassung:
Dein Problem ist die Frequenzauflösung. Wenn du diese verbesserst, triffst du die exakte Frequenz besser und bekommst auch bessere Ergebnisse! Ich empfehle dir also wirklich, diesen "einfachen" Weg zu gehen.



PS. Mein Programm liefert natürlich von der Phase unabhängige Amplituden! ;)

Im übrigen liefert mein Programm auch deine 200 Hz exakt + Phase + Amplitude! Vorausgesetzt natürlich dein Zeitfenster ist nicht länger als dein Signal (0,05 s) siehe Bild ;)

[stoxx]

Hi Peter,

Dein Text trägt zum Verständnis bei, trotzdem funktioniert es nicht :-)

ich hab jetzt mal Deine Normierung unten beim Goertzel Button rausgenommen, so dass man auch eine Frequenz eingeben kann, ohne dass sie verändert wird.

ich gebe ein
300 Hz
3 Amplituede
45 Phase
48000 Abtastfrequenz
0,02 Signaldauer

er findet dann erstmal 288
das änder ich auf 300 und drücke den Goertzel Button

errechnen tut er:

0,0038360000 Amplitude
450 Phase


Hat es einen bestimmten wichtigen Grund, warum Du die den Zähler im Goertzer rückwärts zählst?
(btw wäre besser von Fenstergröße-1 bis null zu iterieren)

[Peter1999]

Das kann ja auch nicht gehen ;)

Wenn du mit 48000 Punkten / Sek eine Signaldauer von 0,02 Sek vorgibst, dann sind das nur 960 Stützstellen! Die Fenstergröße muss dann also kleiner als 960 sein. (Es wundert mich, dass das nicht in eine Exception läuft, aber egal, da müsste ich jetzt mal den Debugger anwerfen). Nimm mal 800 und schon geht's. Sowas habe ich hier in dem Beispiel noch nicht abgefangen.

Das ist ja gerade das Problem bei allen Frequenzbereichsverfahren: höhere Frequenzauflösung benötigt immer längere Signale.