Phasenverschiebung zw. 2 Signalen bestimmen

[moppi (Gast)]

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

[Ronny (Gast)]

Hmmm...Nulldurchgänge angucken? ;)

[Detlef _. (detlef_a)]

>>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

[moppi (Gast)]

Ich habe zum Görtzel folgenden C-Code gefunden:

float goertzel(float freq, int size, const float x[])
{
    int i;
    float coeff;
    float s, s_prev1 = 0.0f, s_prev2 = 0.0f;

    coeff = 2.0f * cosf(2.0f * M_PI * freq);

    for (i = 0; i < size; i++) {
        s = x[i] + (coeff * s_prev1) - s_prev2;
        s_prev2 = s_prev1;
        s_prev1 = s;
    }

    return (s_prev1 * s_prev1) + (s_prev2 * s_prev2)
         - (s_prev1 * s_prev2 * coeff);
}

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?

[Detlef _. (detlef_a)]

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

[Gerrit B. (gbuhe)]

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

[moppi (Gast)]

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

[moppi (Gast)]

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:

realW = 2.0*cos(2.0*pi*k/N);
imagW = sin(2.0*pi*k/N);

d1 = 0.0;
d2 = 0.0;
for (n=0; n<N; ++n)
  {
    y  = x(n) + realW*d1 - d2;
    d2 = d1;
    d1 = y;
  }
resultr = 0.5*realW*d1 - d2;
resulti = imagW*d1;

demnach:

amplitude = sqrt(resultr*resultr+resulti*resulti);
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

[Gerrit B. (gbuhe)]

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

[Jochen W. (moppi)]

Hallo Gerrit, lang nix gelesen!


(Gast) Moppi =/= moppi. Btw.: irgendwer hat sich heut mein "vergessenes" 
Passwort zusenden lassen....
Schon seltsam.



73 de Jochen!

[Karsten Z. (moppi83)]

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!

[Karsten Z. (moppi83)]

Ich bleibe mal beim Goertzel: DA komme ich einfach nicht weiter.
Hier mal mein Delphi-Code:

     amplitude:=StrToInt(Edit1.Text);  //Amplitude abfragen
     scan_freq:=StrToInt(Edit2.Text);  //Scan-Frequenz = Soll-Frequenz
     sample_freq:=StrToInt(Edit3.Text); //Sample-Frequenz
     N:=round(sample_freq/scan_freq);   //Anzahl Messwerte

     SetLength(soll,5*N+1);             //Arraylänge bestimmen

     for i:=0 to 5*N do
     begin
          soll[i]:=amplitude*sin(2*pi*i/N);  //5Sinus-Schwingungen erzeugen
     end;
     

    //goertzel!
     k:=0.5+(N*scan_freq/sample_freq);
     w:=(2*pi/N)*k;
     cosine:=cos(w);
     sine:=sin(w);
     coeff:=2*cosine;
     q1:=0;
     q2:=0;
     for i:=0 to 5*N do
     begin
          q0:=coeff*q1-q2+soll[i];
          q2:=q1;
          q1:=q0;
     end;
     real_value:=q1-q2*cosine;
     imag_value:=q2*sine;
     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