Hallo, ich habe folgendes Problem: Ich habe Messwerte von 2 Sinusschwingungen; Diese sind gleicher Frequenz, jedoch phasenverschoben und mit unterschiedlicher Amplitude; Nun ist die Frage, wie ich die Phasenverschiebung und die Amplitude der Signale rausbekomme. Das Ganze ist zeitkritisch, d.h. die normale DFT fällt schonmal raus. FFT geht durch die Anzahl der Messwerte nicht. Gibt es quasi eine DFT mit nur einer Frequenz (diese stimmt ja überein)??? Grüße Moppi
>>Gibt es quasi eine DFT mit nur einer Frequenz (diese stimmt ja überein)??? Ja, das ist der Goertzel Algorithmus. >>Hmmm...Nulldurchgänge angucken? ;) Schlecht bei schwachem SNR und nicht subsamplegenau. Cheers Detlef
Ich habe zum Görtzel folgenden C-Code gefunden:
1 | float goertzel(float freq, int size, const float x[]) |
2 | {
|
3 | int i; |
4 | float coeff; |
5 | float s, s_prev1 = 0.0f, s_prev2 = 0.0f; |
6 | |
7 | coeff = 2.0f * cosf(2.0f * M_PI * freq); |
8 | |
9 | for (i = 0; i < size; i++) { |
10 | s = x[i] + (coeff * s_prev1) - s_prev2; |
11 | s_prev2 = s_prev1; |
12 | s_prev1 = s; |
13 | }
|
14 | |
15 | return (s_prev1 * s_prev1) + (s_prev2 * s_prev2) |
16 | - (s_prev1 * s_prev2 * coeff); |
17 | }
|
In die Funktion kommt also die Frequenz, die Größe des Arrays und das Array selbst. Ein Float-Wert wird zurückgegeben, welche die Amplitude bei der gegebenen Frequenz darstellt, richtig? Wie komme ich aber zur Phase?
Ja, es scheint irgendwas quadratisches zurückzukommen. Die Phase hat mit dem Verhältnis der s_prev1 und s_prev2 zu tun. Code ist nicht übersetzbar und Zeug wie s_prev2 = s_prev1; s_prev1 = s; ist nicht akzeptabel, such Dir lieber ne andere Implementation. Englische Wikipedia erklärt Goertzel relativ erschöpfend. Cheers Detlef
Hallo Moppi, durch eine Quadraturdemodulation bekommst Du abtastwertgenau die Amplituden- und Phaseninformation. Eines Deiner Sinussignale scheint ja Dein Stimulus zu sein, nennen wir es S_out und dürfte eine sehr gute Qualität haben. Dein in Phase und Amplitude verändertes Detektionssignal nennen wir S_in. Für die Ermittlung der Amplituden- und Phasenänderung von S_in zu S_out benötigen wir noch ein um 90° verschobenes S_out, was ich mit S_out(-90°) bezeichne. Nun folgt die sogenannte Quadraturdemodulation: x = S_out * S_in y = S_out(-90°) * S_in Da bei dieser Demodulation in beiden Pfaden auch ein Signalanteil mit doppelter Frequenz entsteht, ist ein Tiefpaßfilter zur Unterdrückung dieses Anteils nötig. Das kann ganz einfach als IIR-Filter 1.Ordnung implementiert werden (zwei Multiplikationen und eine Addition, bei geschickter Implementierung auch mit Bitshifts statt Multiplikationen realisierbar). Da das Ergebnis nun allerdings kartesisch vorliegt (rechtwinklige Koordinaten x und y), mußt Du dann die Umrechnung in polare Werte(Amplitude und Phase) vornehmen (Hast Du einen Mikrocontroller oder DSP im Einsatz?): Amplitude = wurzel(x^2+y^2) Phase = artan(y/x). Diese Auswertung ist Subsample-genau (nahezu beliebig, wenn bei der Rechnung aufgepaßt wird) und funktioniert auch mit extrem kleinen Signal-Rausch-Abständen. Auf diese Weise arbeiten nämlich Funkempfänger sehr präzise mit sehr kleinen SNRs. Viele Grüße! Gerrit, DL9GFA
Hallo, ich programmiere eine Interfacekarte, die in einem PC hängt (leider mit Delphi, kann ich aber nicht ändern); Ich habe jetzt im Netz gesucht und recht wenig zur Quadraturdemodulation gefunden. Hört sich natürlich erstmal sehr gut an. Danke für den Beitrag Gerrit! Da ich aus der Automatisierungstechnik komme und eigentlich wenig Ahnung von Filtern habe bräuchte ich noch ein paar mehr Infos... . S_in wie auch S_out liegen ja diskreten Abtastwerten vor. Die Rechnung: x = S_out * S_in y = S_out(-90°) * S_in mache ich für jeweils einen Abtastzeitpunkt, richtig? (Das Gerüst drumherum ist mir leider noch nicht richtig ersichtlich). Gibts eine gute URL dafür? Grüße Moppi
Ich verfolge mal auch noch den Goertzel-Ansatz: http://www.mstarlabs.com/dsp/goertzel/goertzel.html Hier ist Goertzel eigentlich recht gut erklärt. Code gibt es auch:
1 | realW = 2.0*cos(2.0*pi*k/N); |
2 | imagW = sin(2.0*pi*k/N); |
3 | |
4 | d1 = 0.0; |
5 | d2 = 0.0; |
6 | for (n=0; n<N; ++n) |
7 | {
|
8 | y = x(n) + realW*d1 - d2; |
9 | d2 = d1; |
10 | d1 = y; |
11 | }
|
12 | resultr = 0.5*realW*d1 - d2; |
13 | resulti = imagW*d1; |
demnach:
1 | amplitude = sqrt(resultr*resultr+resulti*resulti); |
2 | phase = arctan2(resulti, resultr); |
Laut Erklärung " The results are the real and imaginary parts of the DFT transform for frequency k." ist k die Zielfrequenz. Entsprechend bezeichnet n dann die Zahl der Abtastwerte (die Schleife geht ja auch von 0 bis N). Dabei kann man nur eine Periode abtasten, richtig? Funktionieren tuts leider noch nicht.... . Grüße Moppi
Hallo Moppi, am besten kann man es in der numerischen Simulation zeigen und optimieren. Ich verwende gern Scilab dazu (http://www.scilab.org). Anbei findest Du ein Scilab-Skript, das die nötige Mathematik umsetzt, sowie das Ergebnis als Grafik. Ich habe auf S_out einen Amplitudenfaktor von 0.01 und einen Phasenfehler von 30°angewendet, die auch prompt detektiert werden ;o). Mein Skript ist massiv dokumentiert; spiele doch einfach etwas damit herum. Dazu mußt Du Dir nur das freie Scilab installieren und im Dateimanager einen Doppelklick auf moppi01.sce durchführen. Mit <Strg>+<l> oder über Menü des gestarteten Editors kannst Du das Skript ausführen lassen. In der Grafik ist das Einschwingen auf die korrekten Werte zu sehen. Man kann durch Vergrößerung des Filterkoeffizienten die Grenzfrequenz nach oben schieben und damit ein deutlich schnelleres Einschwingen erzielen. Die Unterdrückung der unerwünschten Mischprodukte wird dann aber geringer, was zu größeren Fehlern führt. Das hängt direkt von Deinem Verhältnis Abtastfrequenz zur Trägerfrequenz ab. Mit Filtern höherer Ordnung, die eine hohe Tiefpaß-Grenzfrequenz ermöglichen und trotzdem gleichzeitig eine hohe Unterdrückung des Mischproduktes bei der doppelten Frequenz ermöglichen, kann ein Einschwingen innerhalb weniger Abtastwerte erreicht werden. Zu dieser Optimierung mußt Du Dich dann aber intensiver mit digitalen Filtern auseinander setzen. Ich hoffe das hilft. Viele Grüße! Gerrit, DL9GFA
Hallo Gerrit, lang nix gelesen! (Gast) Moppi =/= moppi. Btw.: irgendwer hat sich heut mein "vergessenes" Passwort zusenden lassen.... Schon seltsam. 73 de Jochen!
Hallo Jochen, sorry, ich dachte ich hätte hier noch nen Account... . Dass es dich mit dem gleichen Nickname gibt wusste ich nicht.... . Grüße Karsten PS: schöner Nickname!
Ich bleibe mal beim Goertzel: DA komme ich einfach nicht weiter. Hier mal mein Delphi-Code:
1 | amplitude:=StrToInt(Edit1.Text); //Amplitude abfragen |
2 | scan_freq:=StrToInt(Edit2.Text); //Scan-Frequenz = Soll-Frequenz |
3 | sample_freq:=StrToInt(Edit3.Text); //Sample-Frequenz |
4 | N:=round(sample_freq/scan_freq); //Anzahl Messwerte |
5 | |
6 | SetLength(soll,5*N+1); //Arraylänge bestimmen |
7 | |
8 | for i:=0 to 5*N do |
9 | begin |
10 | soll[i]:=amplitude*sin(2*pi*i/N); //5Sinus-Schwingungen erzeugen |
11 | end; |
12 | |
13 | |
14 | //goertzel! |
15 | k:=0.5+(N*scan_freq/sample_freq); |
16 | w:=(2*pi/N)*k; |
17 | cosine:=cos(w); |
18 | sine:=sin(w); |
19 | coeff:=2*cosine; |
20 | q1:=0; |
21 | q2:=0; |
22 | for i:=0 to 5*N do |
23 | begin |
24 | q0:=coeff*q1-q2+soll[i]; |
25 | q2:=q1; |
26 | q1:=q0; |
27 | end; |
28 | real_value:=q1-q2*cosine; |
29 | imag_value:=q2*sine; |
30 | measured_amplitude:=sqrt(real_value*real_value+imag_value*imag_value); |
Ist quasi die Implementierung nach http://www.embedded.com/story/OEG20020819S0057 Komischerweise bekomme ich nie die richtige Amplitude angezeigt... . Erkennt jemand einen Fehler? Grüße Karsten
Bitte 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.