Amplitude und Phase einer festen Frequenz bestimmen.

[Mr. Q. (lunza)]

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

[branadic (Gast)]

Hallo,

wie wäre es mit dem Einsatz eines AD8302 (RF / IF Gain Phase Detector )? 
Dann hast du die Problematik mit dem zeitkritisch erschlagen.

branadic

[Toni (Gast)]

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?

[Detlef _. (detlef_a)]

Goertzel 6MHz auf nem uC ist schon tapfer. Gibt bestimmt hochgetaktete 
uCs, mit den mir bekannten geht das allerdings nicht.

Cheers
Detlef

[Mr. Q. (lunza)]

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 (Gast)]

So geht's am schnellsten:

Amplitude mit Spitzenwertdetektion
Phase mit Nulldurchgangsinterpolation

[Mr. Q. (lunza)]

chris schrieb:
> So geht's am schnellsten:
>
> Amplitude mit Spitzenwertdetektion
> Phase mit Nulldurchgangsinterpolation

Wie soll dass funktionieren mit nur 6 Abtastwerten pro periode?

[Mr. Q. (lunza)]

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?

[Strubi (Gast)]

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

[Mr. Q. (lunza)]

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

[Detlef _. (detlef_a)]

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

[Mr. Q. (lunza)]

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

[Detlef _. (detlef_a)]

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

[Mr. Q. (lunza)]

@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

[Mr. Q. (lunza)]

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

[Mr. Q. (lunza)]

ach momentchen, sind ja Matrizen... ich muss das mal eben bearbeiten.

[Mr. Q. (lunza)]

Ich habe nun nur noch das Problem, mit der anzupassenden Verstärkung für 
8,3Mhz. wie berechne ich die?

[Mr. Q. (lunza)]

Noch eine Frage:

Muss das Signal auf 1 Normiert sein? Ansonsten würde ich direkt den 
ADCwert nehmen.

[Mr. Q. (lunza)]

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.

[Fabian H. (Firma: keine) (eimer)]

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

[Detlef _. (detlef_a)]

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