Goertzel Algorithmus (Frequenz + Phase)

[stoxx]

es funktioniert auch mit 800 nicht, warum glaubt mir niemand
(wie Du schon sagst, wenn die Fenstergröße zu groß ist, rechnet er einfach ins leere weiter ..)


ich überlege gerade, ob es entscheidend ist, rückwärts zu zählen.
aus einem Buch über Numerische Mathematik bei google books hab ich entnommen, das es wohl so ist.
Da drin sieht man auch, ich nehme an, dass der Autor Ahnung hat, dass man um den Real Teil auszurechnen, cosine mit Q1 multiplizieren muss. in dem Link ist es gerade vertauscht (und falsch) und führt also auch zu diesen Problemen.

http://www.eetimes.com/design/embedd...tzel-Algorithm

Trotzdem geht aber die Phase immer noch nicht

Vielleicht probier ich das mal, mit dem Index rückwärts laufen lassen.

In dem Buch für Numerische Mathematik steht auch, dass der Algorithmus nicht gutartig wäre, was auch immer das heißt.
Es wird der Algorithmus von Reinsch empfohlen.
der macht bei kleinem Epsilon ? weniger Fehler .. ob das wichtig ist?
oder sind das trotzdem nur kleine Fehler?
was denkst Du?

[Peter1999]

Doch, ich glaube dir schon .
Man macht sich ja nicht freiwillig das Leben schwer .

Ich kann das Problem nur nicht nachvollziehen, da es ja mit deinen Werten funktioniert, wenn ich mit einer Fenstergröße von 800 Punkten rechne (siehe Bild).

Ob man rückwärts zählen muss (oder nur kann) weiß ich nicht, ich habe den Algorithmus nicht erfunden. Und ich glaube, es ist jetzt auch zu spät, um noch schnell die Herleitung zu verstehen (Habe auch gerade kein vollständiges Lehrbuch zur Hand ).

Willst du nicht einfach die Variante aus meinem Beispiel übernehmen und gut ist's?
Ich habe jedenfalls noch keinen Fall gefunden (der plausibel ist) der nicht funktioniert.

[stoxx]

tut mir Leid .. ich hatte mich vorhin bei den Zahlen verschrieben.

(wegen unzureichend getesteter Software wird irgendwann sicher nochmal die Welt untergehen

gib mal ein
Frequenz 300
Amplitude 3
Phase 50
abtastfrequenz 8000
Signaldauer 0,04 oder so

dann Du sehen

[Peter1999]

Das geht auch !!

8000 Punkte / Sek * 0,04 Sek = 320 Punkte

--> gewählt Fenstergröße 300 Punkte (weil < 320)

--> gefundene Frequenz 292,333 Hz, das ist natürlich nicht so gut (daher mein Vorschlag mit der Optimierung)

--> Häckchen entfernt: 300 Hz vorgegeben (die würde die Optimierung dann ja auch automatisch finden)

und schwups: Korrektes Ergebnis für Ampl + Phase!

[stoxx]

schick mal bitte einen Screenshot ..
oder hat Intel jetzt wieder einen Rechenfehler in seiner CPU?
Da der Algo noch "falscher" ist bei Dir, kanns eigentlich gar nich gehen

[stoxx]

so .. diesen Fehler hab ich jetzt gefunden .. interessant.
Dein altes Delphi7 compiliert (zum Glück) anders .. sonst hätte man das Problem wohl nicht gleich gefunden .. bzw lag es genau daran, was im Buch besprochen wurde, dass es kein "guter" Algorithmus wäre, wegen Fehleranfällgikeit. (hier Fehlerfortpflanzung)

soo .. Deine Normierung werde ich mir also jetzt mal genauer ansehen !

hab vielen Dank ! .. meine Nacht ist gerettet
(Die klammer in der For Schleife)


Mit einer Klammer, die Du sicherheitshalber auch noch setzen kannst, geht es auch jetzt bei mir. erstaunlich..

Delphi-Quellcode:

  q1 := 0;
  q2 := 0;
  for i := Fenstergroesse downto 0 do  
  begin                                
    q3 := Signal[i]+ ( alpha*q1-q2 );
    q2 := q1;
    q1 := q3;
  end;


  re := q1-(alpha*q2)/2;
  imag := q2 * sin(2 * Pi * Frequenz);

  re := re / Fenstergroesse * 2;
  imag := imag / Fenstergroesse * 2;

im Prinzip kann man sich aber bei diesem Algo die Fenstergröße sparen ...

[stoxx]

Hi Peter,

darf ich Dich noch bissl nerven und Dir nich ein paar Fragen stellen?

wenn ich eine Frequenz von 300 Hz erzeuge, und dann eine Abfrage erstelle, dann kommt für Frequenzen von 750 bis 850 folgendes raus:
Was mich dabei stört, der der Amplitudenwert schwankt zwischen 0,82 und 5 .. obwohl da gar keine Frequenz vorhanden ist.
Meine Vermutung ist jetzt, dass es mit der unpassenden Fenstergröße zusammenhängt.
Was ist also mathematisch die korrekte nächst kleinere Fenstergröße um einen minimalen Fehler zu bekommen?
Hat das was mit diesem Fehler zu tun, von dem hier gesprochen wird?


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

Und dann noch eine Frage ..
Und weißt Du, ob man den Algorithmus umbauen könnte, die Amplitude zu bestimmen, wenn ich eine Frequenz auf einer bestimmten Phase suche?
Wir haben den Fall, dass wenn im Signal 300 Hz drin sind, mit einer stark abweichenden Phase, dass man dann, wenn man exakt diese Frequenz abfragt, nicht mehr das korrekte Amplitudenmaximum errechnet werden kann. scheinbar fragt der Algo immer Phase Null ab..

Danke Dir schon mal wieder

Delphi-Quellcode:

frequenz: 756 ;real: -2,91582615873547 ;imag: -2,57019064576388 ;Amp: 3,89
frequenz: 760 ;real: -2,33132644235365 ;imag: -4,34704698053365 ;Amp: 4,94
frequenz: 764 ;real: -0,578012729919608 ;imag: -4,99372967714084 ;Amp: 5,03
frequenz: 768 ;real: 1,03 ;imag: -4,03398585057834 ;Amp: 4,17
frequenz: 772 ;real: 1,31 ;imag: -2,18667224046469 ;Amp: 2,55
frequenz: 776 ;real: 0,07 ;imag: -0,815913068896486 ;Amp: 0,82
frequenz: 780 ;real: -1,74999260357873 ;imag: -0,912931983818342 ;Amp: 1,98
frequenz: 784 ;real: -2,79015995586682 ;imag: -2,36906710716151 ;Amp: 3,67
frequenz: 788 ;real: -2,29737731999792 ;imag: -4,07268794135492 ;Amp: 4,68
frequenz: 792 ;real: -0,65868231160542 ;imag: -4,74592185313099 ;Amp: 4,8
frequenz: 796 ;real: 0,9 ;imag: -3,89177044592188 ;Amp: 4
frequenz: 800 ;real: 1,23 ;imag: -2,150048665788 ;Amp: 2,48
frequenz: 804 ;real: 0,1 ;imag: -0,807807734222314 ;Amp: 0,82
frequenz: 808 ;real: -1,63655987839042 ;imag: -0,838429931342038 ;Amp: 1,84
frequenz: 812 ;real: -2,67535127160256 ;imag: -2,18643363723128 ;Amp: 3,46
frequenz: 816 ;real: -2,26671759674278 ;imag: -3,82308911553845 ;Amp: 4,45
frequenz: 820 ;real: -0,733016523832668 ;imag: -4,52076309233921 ;Amp: 4,58
frequenz: 824 ;real: 0,78 ;imag: -3,76345212910821 ;Amp: 3,85
frequenz: 828 ;real: 1,15 ;imag: -2,11850443555382 ;Amp: 2,41
frequenz: 832 ;real: 0,12 ;imag: -0,802207476643857 ;Amp: 0,82
frequenz: 836 ;real: -1,53194381311629 ;imag: -0,771490400560849 ;Amp: 1,72
frequenz: 840 ;real: -2,56954942411791 ;imag: -2,01972832187638 ;Amp: 3,27
frequenz: 844 ;real: -2,23856498360031 ;imag: -3,59452992539954 ;Amp: 4,24
frequenz: 848 ;real: -0,801801997439787 ;imag: -4,31466545944052 ;Amp: 4,39

[Peter1999]

Hallo,

ich versuchs nochmal ganz von Anfang .
Ich muss aber dazu sagen, ich bin KEIN Mathematiker, es kann also alles falsch sein, was ich hier von mir gebe.

1. Du brauchst unbedingt eine Fensterung! Ich weiß ja nicht, wie lang deine Signale sind, aber die DFT (bzw.) der Gortzel benötigt N^2 Rechenoperationen. Wenn du da das ganze Signal als Fenster nimmst, bist du bei 5min*60sek/min*48000samples/sek ganz schnell am Ende des machbaren. HINWEIS: Auch das ganze Signal ist ein Zeitfenster, eben ein Rechteckfenster der max. Größe.
Daher wurde ja die FFT erfunden mit "nur" NlogN Rechenoperationen.

2. Der Goertzel ist KEIN Verfahren zur Frequenzbestimmung! Es ist ein einfacher Algorithmus zur Transformation aus dem Zeit- in den Frequenzbereich. Was hier gemacht wird ist eigentlich eine "Vergewaltigung" des Verfahrens. ANNAHME: es existiert genau EIN (oder mehrere, dann aber im Frequnzbereich weit von einander eintfernte, typischer Fall: Analogtechnik) Signal, einer Frequenz, dessen Amplitude deutlich größer ist als die der Störungen. Was wir hier machen ist also: Suche den gesamten Amplitudengang ab bis du das Maximum gefunden hast und merke dir die zugehörige Frequenz.

3. Da du es IMMER mit einem Zeitfenster zu tun hast (größe hin oder her), hast du es auch IMMER mit den typischen Problemen zu tun. Das wäre in erster Linie die Periodisierung aufgrund der Disketisierung. Daraus folgt der "Abschneidefehler", der mit verschiedenen Fensterfunktionen reduziert werden kann. Häufigstes Beispiel: Hanning. Du hast es also IMMER mit Seitenbändern zu tun. (Schau dir dazu mal den Amplitudengang eines normierten Rechteckfensters an (numerisch oder Laplacetransformierte...ist egal), dann siehst du das Abklingen mit dem Abstand von der "echten" Frequenz.

4. Ich vermute, was du hier geplottet hast ist genau der Effekt der Nebenbänder. Das kann ich aber nicht sicher sagen, da ich deine Parameter nicht genau kenne. Ich habe in meinem Programm mal etwas gespielt, da lande ich bei höheren Frequenzen in der Größenordnung 1% der "echten" Amplitude. Vielleicht ist bei dir aber auch noch ein Fehler im Quelltext

5. Nein, es ist nicht möglich, zu einer Phase die Amplitude zu finden. Jedenfalls ist es nicht sinnvoll. Dann müsstest du den gesamten Phasengang bestimmen (so wie wir es jetzt für die Amplitude machen), bei einer bestimmten Phasenüberschreitung anhalten und dir dann die Amplitude und Frequenz ausgeben lassen. Wenn man so etwas machen will, ist man wieder bei der FFT, die da deutlich schneller ist.

6. Den Goeztzel würde ich eigentlich nur verwenden, wenn ich die Frequenz schon vorher kenne (z.B. Sender-Frequenz) oder wenn ich zwischen zwei Spektrallinien einer FFT eine Interpolation durchführen möchte. Dann kann man sich schnell und mit wenig Aufwand ein paar Stützstellen berechnen, aber KEIN ganzes Spektrum!

7. Du beachtest aber schon die ganze Zeit deinen max. Frequenzinhalt? Es gilt natürlich Shannon! Also Abtastfrequenz / 2 ist deine maximal auswertbare Frequenz. Danach greift wieder die Periodisierung des Spektrums und du bekommst die berühmten Spiegelfrequenzen. Also niedrige Frequenzen tauchen bei höheren Frequenten wieder auf. (vielleicht ist das auch das Problem deiner Amplituden)

8. Die Frequenzauflösung ist immer fest mit der Fenstergröße verbunden! Wenn du also nicht genau auf der richtigen Frequenz landest, geht deine Amplitude nach unten und die Phase läuft dir weg. delta f = 1 / T .

9. Willst du auch diesen Fall absichern, musst du, wie gesagt, ganz anders vorgehen und aus dem Verlauf des Amplitudengangs das Maximum und dem Phasengang die zugehörige Phase bestimmen. (Damit wären wir bei der Eigenfrequenzbestimmung und der modalen Dämpfung ) Das ist dann aber ein ganz anderes Kapitel.

10. Probier mal, ob dein Problem mit den Amplituden auch in meinem Programm auftritt (wenn du die Sachen wie Shannon & Fenstergröße beachtest).

11. Ich habe mein Programm mit Delphi 2009 UND Delphi 7 geschrieben und kompiliert, es funktionieren beide EXEn gleich.

12. Was mir heute morgen eingefallen ist: Es ist natürlich logisch, das die Phasen falsch werden, wenn man am falschen Ende des Zeitfensters anfängt. Das ist dann ja von der Signalfrequenz und der Fenstergröße abhängig. Und deine Buchseite sagte ja auch, man solle rückwärts Summieren.


So, damit habe nun eigentlich alles aufgeschrieben, was ich dir so aus dem Kopf sagen kann. Viel mehr kann ich dir dann sicher auch nicht helfen. Ich kann dir aber noch ein altes DDR-Buch empfehlen: "Experimentelle Festkorpermechanik", VEB Fachbuchverlag, Leipzig 1986. Da ist die Frequenzbereichstransformation schön übersichtlich mit vielen Bildchen für Dummies (damit meine ich mich) erklärt

Viel Glück

[stoxx]

ja .. es ist eine Vergewaltigung des Verfahrens.
ABER
ich kenn mich mit FFT nicht aus, und ich kann nur nachplappern, was ich gelesen habe.
Man braucht bei der FFT immer eine Fensterung .. Genau die entfällt hier.
(Vorteil, da ich einen ganzen Bereich analysieren kann)

Deswegen gibr es auch kein "abflachen" am Fensterrand.
Außerdem ist man bei FFT auf Frequenzen mit Potenzen von 2 beschränkt.
Stimmt das eigentlich, oder wie umgeht man das dann?

Zeitkritisch sind die Sachen erstmal nicht, von daher sind die Rechenoperationen egal.

> 1% der "echten" Amplitude

ja .. das ist zuviel. Ein Algorithmus sollte nicht abhängig von dem Verhälltnis Samplerate/Targetfrequenz sein.

ich werde jetzt folgendes probieren, ich werde die nächst kleinere Fenstergröße wählen, einmal links ansetzen, einmal rechts ansetzen, und den Mittelwert davon bilden.
mal sehen, was dabei rauskommt, somit habe ich alle daten aus dem Fenster, aber mit besserer Amplitude.

> Dann müsstest du den gesamten Phasengang bestimmen

das heißt? .. wüsstest Du denn, wie das geht? würde das auch mit Goertzel gehen?


> aber KEIN ganzes Spektrum!

wie gesagt, da ich den Algo zweckentfremden möchte, ich aber keine Ahnung von FFT habe, würde ich sie nie anwenden, wenn ich nicht das Teil komplett verstanden, also selbst entwickelt habe.
Von daher fällt FFT flach, sonst kommt nur Murcks dabei raus.

Ich hab z.b. jetzt ein Array mit daten, und wenn ich die ersten 2 Elemente weglasse, kann ich über den Phasenwinkel die Verschiebung berechnen, und komme auf "2 Elemente verschoben" .. das ist toll und macht mich glücklich

oder andere Frage .. könnte denn die FFT das Gleiche?
Also zu JEDER Frequenz die Amplitude und Phase und zu bestimmen?
In EINEM Rechteck.

Danke für Deine ausführlichen Erklärungen

[TiGü]

Zitat von stoxx:

ja .. es ist eine Vergewaltigung des Verfahrens.
ABER
ich kenn mich mit FFT nicht aus, und ich kann nur nachplappern, was ich gelesen habe.
Man braucht bei der FFT immer eine Fensterung .. Genau die entfällt hier.
(Vorteil, da ich einen ganzen Bereich analysieren kann)

Du brauchst auch eine Fensterung. Zur Zeit hast du ein riesiges Rechteckfenster. Begründung siehe 1. von Peter1999.

Zitat von stoxx:

Deswegen gibr es auch kein "abflachen" am Fensterrand.

Unkritisch, denn so wie ich verstanden habe, hast du eine handvoll Frequenzen in einer überschaubaren Bandbreite.
Ich rate dringend sich mit elementarer digitaler Signalverarbeitung auseinanderzusetzen. Findet sich genug im Netz dazu.

Zitat von stoxx:

Außerdem ist man bei FFT auf Frequenzen mit Potenzen von 2 beschränkt.
Stimmt das eigentlich, oder wie umgeht man das dann?

Das stimmt so nicht!
Der Algorithmus arbeitet effizient wenn die Anzahl der Stützstellen / Abtastwerte eine Potenz von zwei ist.

Zitat von stoxx:

Ein Algorithmus sollte nicht abhängig von dem Verhälltnis Samplerate/Targetfrequenz sein.

Mein Nachrichtentechnik-Studium liegt schon ein bissel zurück, aber was du forderst sehe ich nicht.
Vielleicht missverständlich ausgedrückt?
Kann ja schließlich kein 20 MHz Signal mit 500 Hz abtasten!
Liest sich so, als ob ein Auto unabhängig von Sprit immer gleich gut, weit und schnell fahren soll.

Zitat von stoxx:

wie gesagt, da ich den Algo zweckentfremden möchte, ich aber keine Ahnung von FFT habe, würde ich sie nie anwenden, wenn ich nicht das Teil komplett verstanden, also selbst entwickelt habe.
Von daher fällt FFT flach, sonst kommt nur Murcks dabei raus.

Schau dir das mal an:
http://www.delphipraxis.net/87967-ff...gorithmus.html

Ansonsten Google -> Delphi + FFT

Des Weiteren gibt es genügend Beispiele für Implementationen in anderen Programmiersprachen.

Wozu brauchst du das eigentlich?