Hallo,
ich möchte die Amplitude und Phase einer festen Frequenz bestimmen.
Das Signal hat eine Frequenz von 1Mhz und wird mit 6 Mhz abgetastet. Ich
gehe davon aus, das durch eine angemessene analoge filterung keine
Oberwellen o.ä mehr vorhanden sind.
Das Referenzsignal wird nicht abgetastet, die Phase ist jedoch bekannt.
Nun möchte ich die Amplitude des Messignal bestimmen und die
Phasenverschiebung zu dem anderen Signal.
Das Problem ist dass ganze wird sehr zeitkritisch in nem uC stattfinden.
Nun ist die überlegung über einen Goertzel Filter Real und Imaginärteil
zu bestimmen und in Amplidude und Phase zurück zu rechnen.
Da ich aber wirklich sagen kann das die Frequenz immer stabil ist, gäbe
es da noch andere möglichkeiten? Je schneller desto besser.
Gruß
Mr. QB
Der Goertzel ist schonmal eine gute wahl, wenn Du nur eine Frequenz
haben möchtest. In zeitkritischen Systemen ist er (abhängig vom DSP)
recht schnell verarbeitet, einzig die Berechnung von Betrag und Phase
könnte von der Zeit her etwas kniffelig werden, schau dir dazu mal den
Cordic - Algorithmus an.
Welchen DSP verwendest Du?
Hallo,
Erstmal danke für die schnellen Antworten. Also den ad8303 hatte ich mir
vorher schonmal angesehen, ich bin mir aber sehr unsicher ob mir die
Auflösung in Betrag und Phase recht, da ich mich sehr weit im
Grenzbereich befinde.
Den Cordic werde ich mir morgen mal ansehen, ob ich an der Stelle
effektiver bin. Einsetzen werde ich nen uC von Texas mit 150Mhz und
fließkommaeinheit.
Der goerzel wird ja auf 1mhz angewendet. Also zu schaffen ist das
allemal, nur fragt sich wie schnell. Ich bin ja für jede Idee offen.
chris schrieb:> So geht's am schnellsten:>> Amplitude mit Spitzenwertdetektion> Phase mit Nulldurchgangsinterpolation
Wie soll dass funktionieren mit nur 6 Abtastwerten pro periode?
Das Problem ist ich benötige einen kontinuierlichen Datenstrom. Meint
ihr ich könnte eine Art IIR Filter verwenden? Oder ist es Außreichend
den Goerzel immer nur über 6 Sample laufen zu lassen?
Hi,
Der Goertzel ist tatsächlich ein IIR-Filter. Ansich nur eine
Nachmodellierung eines angeregten Oszillators. Du kannst nach 6 Samples
den Akkumulator immer wieder löschen und mit den nächsten 6 Samples
wieder los integrieren oder die Koeffizienten so wählen, dass du
faktisch eine Dämpfung hast, und dazwischen den Akku nicht nullen musst.
Weiss halt nicht, womit du dein Filter modellierst , aber es lässt sich
auch von Hand herleiten.
Apropos, wenn die Frequenz gut eingrenzbar ist und keine anderen
Artefaktfrequenzen auftreten (Lowpass!) gibt es auch Tricks mit
Unterabtastung. Dann detektierst du nur die Schwebung mit Goertzel und
rechnest auf die eigentlich gemessene Frequenz zurück.
Gruss,
- Strubi
Mahlzeit,
Welchen Einfluss haben denn die Koeffizienten? Also ich halte mich da an
die Mathemathik und berechne folgendermaßen:
_u32 k;
_f32 w;
k = (_u32)(0.5+((N*TARGET_FREQ)/SAMPLE_RATE));
w = ((2*pi)/N)*k;
cosine = cos(w);
sine = sin(w);
coeff = 2*cosine;
Was meinst du mit nicht nullen? Quasi ein gleitender Wert? Das wäre im
Prinzip das Optimalste.
Also die Frequenz ist sehr gut eingrenkbar. Ich habe vor in einer
calibrierung die exakte Frequenz zu erfassen und praktisch auf dies
einzurasten.
Wenn ich ich unterabtaste verschiebe ich ja quasi auch nach "unten", bin
ich dann im endeffekt nicht langsamer?
Wie meinst du das mit Schwebung? der Begriff ist mir neu, wenn es so
effektiver ist wäre es einen Versuch wert.
Gruß Mr.QB
Du kannst die Phase/Amplitude des Signals fortlaufend schätzen, also neu
für jedes reinkommende sample. Das geht mittels eines Kalman Filters.
Das hatte ich dort
Beitrag "Sinus mit 2f < Abtastrate < 4*f: Schwankung herausrechnen"
schon mal breit dargelegt. Angehängt der Matlab Code für Deine Frequenz,
6 Abtastwerte/Welle und der Plot eines Ergebnisses. Mit der Wahl der K
bestimms Du, wie 'nervös' Dein Schätzer ist, also schnell/ungenau oder
langsam/genau. Der Schätzer liefert Dir für jedes sample eine neue Phase
des Systems, die sollte immer um 2*pi/6 weiterspringen. Falls die um
einen anderen Wert weiterrückt kannst Du daraus auch Deine real
vorhandene Frequnz bestimmen.
Pro sample benötigst Du bei dem Vorgehen ca. 7-8Multiplikationen, das
sollte Deine Hardware bei 6Ms/s schaffen.
So würde ich das machen
math rulez!
Cheers
Detlef
PS
>>Das wäre im Prinzip das Optimalste.
'Optimal' läßt sich nicht steigern, das ist ja schon am Optimalsten
clear
%alpha=pi/10;
% Das ist der Drehwinkel des Zeigers
%6 samples pro Welle
alpha=2*pi/6;
n=300;
%ungedämpfter Schwinger
Ap=(1-0*1e-2);
%die Rotations Matrix
A=Ap*[ cos(alpha) -sin(alpha) ; ...
sin(alpha) cos(alpha) ];
%Ini Zustandsgrößen
x=zeros(2,n);
x(1,1)=1;
for(k=2:length(x)) x(:,k)=A*x(:,k-1);end;
% Kalman Systemmodell
C=[0 1];
y=C*x;
%Meßwerte bißchen verrauschen
%y=y+0.1*randn(1,length(y));
A=Ap*[ cos(alpha) -sin(alpha) ; ...
sin(alpha) cos(alpha) ];
P=100*eye(2);
Q=1e1*eye(2);
R=1e3;
% Diagnose Variablen zum Filter
interim=[];
xxd=[];
xd=[0 ; 0];
KK=[];
yyd=[];
% Das sind die stationären Kalman Verstärkungen
ks=[ -0.0533 ; 0.1214];
for k = 1:n,
xxd =[xxd xd];
%PS = A*P*A'+Q;
%K = PS*C'*inv(C*PS*C'+R);
K=ks;
%P = (eye(2)-K*C);
%P = P*PS*P'+K*R*K';
%KK = [KK K];
xds=A*xd;
xd = xds+K*(y(k) -C*xds);
yyd=[yyd C*xd];
%interim=[interim sum(diag(P))];
interim=[interim sum(abs(K))];
end;
%plot([xxd(1,:); x(1,:)].')
%plot([C*xxd ; abs(xxd(1,:)+j*xxd(2,:));y].')
plot(1:n,y,'b.-',1:n,yyd,'r.-')
%plot([abs(x(1,:)+j*(x(2,:)));abs(xxd(1,:)+j*(xxd(2,:)))].')
%plot(interim);
return
Also die Lösung mit dem Kalman gefällt mir sehr.
Die Amplitude und die Phase meines Signals wird aber immer um einen Wert
schwanken den ich messen will, ich hoffe die Genauigkeit reicht,
ansonsten werde ich nach ein Paar mehr samples eine genaue
Zwischenberechnung durchführen und dann einfach immer mal abgleichen.
Ich muss mich noch mit der Mathematik dahinter befassen, aber mir sind
ein paar Dinge noch nicht einleuchtend.
Warum berechnest du
A=Ap*[ cos(alpha) -sin(alpha) ; ...
sin(alpha) cos(alpha) ];
zwei mal?
Deie Berechnung befindet sich ja in:
xxd =[xxd xd];
%PS = A*P*A'+Q;
%K = PS*C'*inv(C*PS*C'+R);
K=ks;
%P = (eye(2)-K*C);
%P = P*PS*P'+K*R*K';
%KK = [KK K];
xds=A*xd;
xd = xds+K*(y(k) -C*xds);
yyd=[yyd C*xd];
%interim=[interim sum(diag(P))];
interim=[interim sum(abs(K))];
y wäre dann mein aktueller Abtastwert oder?
was macht interim?
Wie komme ich auf Amplitude und Phase?
Gruß QB
>>Warum berechnest du ... zwei mal?
Ist egal, kannst Du auch einmal machen
>>y wäre dann mein aktueller Abtastwert oder?
genau
>>was macht interim?
Das steht unter Diagnose Variablen. Da habe ich mir für Tests den
Verlauf der K mitgeloggt. Das ist für die Funktion nicht wesentlich.
>>Wie komme ich auf Amplitude und Phase?
Amplitude: sqrt(xd(1)^2+xd(2)^2) = abs(xd)
Phase: atan2(xd(2),xd(1)); = winkel(xd)
xd ist ein System zweiter Ordnung, ein Vektor mit zwei Komponenten,
Realteil und Imaginärteil eines (komplexen) Vektors. Der Vektor dreht
sich durch xd=A*xd jeweils um 2*pi/6 rad weiter, also nach 6
Abtastwerten ist er einmal 'rum'. Der Ausgangsgleichung dieses Systems
ist yd=C*xd= [0 1]*xd, du siehst also nur eine Komponente dieses
Vektors, die sieht aus wie ein Sinus. Mit den Zeilen xds=A*xd; xd =
xds+K*(y(k) -C*xds);
'schätzt' Du dieses System, Du bestimmst anhand Deiner Meßwerte y die
Zstandsgrößen xd. Erster Schritt: Wie müßte das alte System sich
weitergedreht haben: xds=A*xd; Dann Korrektur dieser 'a priori'
Schätzung mit dem Meßwert, die zweite Gleichung.
Wenn Du die K's festlegst brauchst Du nur diese beiden Gleichungen
auszurechnen.
Cheers
Detlef
@Detlef
vielen dank für die ausführliche Beschreibung, ich bin schon ganz
kribbelig es auszuprobieren.
Ich werde montag aller wahrscheinlichkeit dazu kommen das ganze mal in
c- Code umzusetzen, werde den hier dann mal posten. Bis dahin noch nen
schönen Tag.
Gruß QB
Mahlzeit,
so ist leider ne weile her, das ich an dem Thema arbeiten konnte.
Folgedes:
Ich habe meine Abtastfrequenz dezeit auf 8,3Mhz angepasst. Wenn ich nun
mit deiner Formel eine neue Verstärkung ausrechene:
%PS = A*P*A'+Q;
%K = PS*C'*inv(C*PS*C'+R);
fällt mir auf, das C = {0,1} ist. Somit ist ein Wert des Feldes immer 0.
Kann das richtig sein?
Die Formeln als Quellcode:
xds[0]=KAL_a*xd[0];
xds[1]=KAL_a*xd[1];
xd[0] = xds[0]+KAL_ks[0]*(y -c[0]*xds[0]);
xd[1] = xds[1]+KAL_ks[1]*(y -c[1]*xds[1]);
Auch hier Spielt der hintere teil der Klammer bei C[0] = 0 im Grunde ja
keine Rolle.
Sehe ich das richtig, oder ist da noch der Wurm drin?
Welche Rolle spielt, ob das Messignal eine Offsett hat, oder sollte ich
den vorher Rausrechnen?
Gruß
Mr. QB
Also, mein Filter Läuft, aber eines hast du wohl falsch verstanden oder
ich mich unklar ausgedrückt.
Ich suche nicht die aktuelle Amplidude und Phase eines Messeweres,
sonder einer gesammten Sinusschwingung.
Hi Leute,
ich möchte diesen Thread noch einmal aufgreifen, ich möchte das ganze
eventuell implementieren, verstehe aber ein paar Dinge nicht.
1) Welche Frequenz hat die Schwingung?
[c]alpha=2*pi/6;
A=Ap*[ cos(alpha) -sin(alpha) ; ...
sin(alpha) cos(alpha) ];[\c]
Mit diesem Code wird ein Drehzeiger erzeugt, der zwischen jedem
"Abtastwert" um 120° weiter dreht. Wenn ich das richtig verstehe bildet
der Winkel mit der entsprechenden Abtastzeit die Frequenz.
Bsp:
120° ~ 1ms >> damit ist eine Periode 360°/60°*1ms=6ms lang also 1/6kHz.
Richtig?
2) Wie wird dem Filter mitgeteilt auf welcher Frequenz er schwingen
soll?
ks fällt ja vom Himmel, es ist weder eine Zeit noch eine Frequenz
angegeben. Es werden nur Abtastpunkte erzeugt, ich könnte aber jetzt die
gleichen Abtastwerte für jede beliebige Frequenz erhalten, und einfach
die Zeitbasis ändern und damit eine höhere Frequenz erhalten. Nur unter
der Annahme das die Zeitachse gleich bleibt müsste sich eine höhere
Phasenänderung ergeben oder?
Daraus folgt Frage 3)
Ich möchte die Phasendifferenz von einem Abgetasteten Signal und einem
Referenzsignal bestimmen. Die Eigenfrequenz müsste ich ja
dementsprechend mit dem Kalmanfilter definieren oder?
Ich hoffe der Thread lebt noch :)
Danke
Yo, lebt noch.
>>Richtig?
Richtig, 1/6 kHz. 1ms zwischen zwei Abtastwerten, 6 Werte pro voller
Welle, also Signalfrequenz 1/6 kHz
>>Wie wird dem Filter mitgeteilt auf welcher Frequenz er schwingen
soll?
Durch die A Matrix, 2*pi/6 ist die Freqenz.
>>ks fällt ja vom Himmel, es ist weder eine Zeit noch eine Frequenz
angegeben.
ks hat mit der Frequenz überhaupt nichts zu tun, das sind die
stationären Kalman-Verstärkungen. Das ist der stationäre Wert der 'K'
für die angenommenen Meßrausch- und Eingangsrauschverhältnisse. Das ks
bestimmt wie 'nervös' und schnell oder 'ruhig' und langsam die Schätzung
läuft.
>>für jede beliebige Frequenz erhalten,
yo, das ist bei digitaler Signalverarbeitung immer so. Es ist nur das
Verhältnis von Abtastf. zu Signalf. relevant.
>>Die Eigenfrequenz müsste ich ja
dementsprechend mit dem Kalmanfilter definieren oder?
ja, alpha ist 2*pi*fsignal/fabtast, in obigem Fall also
2*pi*(1/6 kHz)/1kHz=2*pi/6
Cheers
Detlef
Die Angabe einer E-Mail-Adresse ist freiwillig. Wenn Sie automatisch per E-Mail über Antworten auf Ihren Beitrag informiert werden möchten, melden Sie sich bitte an.
Wichtige Regeln - erst lesen, dann posten!
Groß- und Kleinschreibung verwenden
Längeren Sourcecode nicht im Text einfügen, sondern als Dateianhang