Hi,
Mein Problem: Wenn ich einen Sinus mit einer Abtastrate kleiner 4* der
Frequenz des Signals messe, bekomme ich ein Ergebnis, wie es der
Screenshot anzeigt: die gemessene Amplitude schwankt regelmäßig zwischen
der tatsächlichen Amplitude und Null.
Woher das Problem kommt ist mir klar, bei zu wenigen Messpunkten messe
ich das Signal mal wenn es die volle Amplitude hat, bei der nächsten
Messung dann bei etwas weniger usw. bis der Messpunkt das Signal quasi
durchlaufen hat und wieder bei der vollen Amplitude angekommen ist.
Das Problem bei der Sache ist, dass ich das (erste und zweite) Integral
des Signals in Labview berechnen will, schwankende Messwerte für die
Amplitude ergeben dabei eine niedrigere gemessene Amplitude als real
vorhanden ist. Gibt es da eine Methode, um diese Schwankung
herauszurechnen? Sie scheint sich ja als (auf das Signal
multiplizierter) Sinus darzustellen, wenn es also eine Formel gibt mit
der ich aus Signalfrequenz und Abtastrate diese Schwankung ausrechnen
kann sollte das möglich sein, nur habe ich beim Googeln nichts
dergleichen gefunden, deshalb hoffe ich dass mir hier geholfen werden
kann.
Hallo Syndic,
die eigentliche Lösung für Dein Problem heißt digitales Upsampling. Das
ist nicht schwierig, aber erfordert ein bißchen Ahnung von digitaler
Signalverarbeitung. Oftmals ist dies aber nicht erforderlich. Dazu müßte
man aber nochmal genau wissen, was Du machen willst, denn unter den
Begriffen 1./2. Integral kann ich mir nichts vorstellen. Vielleicht:
Etwas bei dem Problem versteh ich nicht so recht: Mehr als zwei, aber
weniger als vier Samples pro Periode sollen bei einer reinen
Sinusschwingung nicht ausreichen um sie zu rekonstruieren? Glückwunsch,
Du hast das Abtasttheorem widerlegt. ;-)
Die Grafik zeigt das abgetastete Signal, ja? Also im Prinzip sollte es
keine kontinuierliche, "durchverbundene" Kurve sein sondern nur einzelne
Punkte, seh ich das richtig? Auffallend ist dass immer zwischen zwei
Punkten das Vorzeichen wechselt. Das würde ja nur gehen, wenn jeweils
ein Sample aus der positiven Halbwelle des Sinus stammt und ein Sample
aus der negativen Halbwelle. Ansonsten müsstest Du auch mal zwei
positive oder zwei negative Abtastwerte hintereinander haben, wie man an
dieser Wikipedia-Grafik sehen kann:
http://upload.wikimedia.org/wikipedia/commons/9/9c/Nyquist_Aliasing.svg
Wenn dem nicht so ist, dann kann Deine Aussage mit Abtastrate >
2*Signalfrequenz wohl nicht stimmen.
Andere Möglichkeit: Das verwendete analoge Tiefpassfilter in der
Signalstrecke hat keine besonders hohe Flankensteilheit, und man müsste
die Abtastrate entsprechend erhöhen (Überabtastung) um das zu
kompensieren. Die Grafik scheint dafür zu sprechen, dass man sehr knapp
an f_Abtast > 2 * f_Signal dran ist (denn f_Abtast = 3 * f_Signal oder
f_Abtast = 3,5 * f_Signal würde definitiv anders aussehen, jedenfalls
wenn alle Komponenten auf der Signalverarbeitungssrecke entsprechend
mitspielen).
Siehe auch http://de.wikipedia.org/wiki/Abtasttheorem#.C3.9Cberabtastung
Das Beispielbild oben ist mit einer Signalfrequenz von 2450 Hz und einer
Abtastrate von 4994,29 sps entstanden, die Abtastrate ist also knapp
mehr als 2*f ;)
markbrandis: das 6. Diagramm in dem Bild das du zitiert hast zeigt recht
gut, was passiert - man hat zwar knapp mehr als 2 Messpunkte pro
Periode, wenn man aber nur die Messpunkte verbindet, bekommt man ein
Bild ähnlich dem, das ich gepostet habe. Ich sage nicht dass ich das
Signal nicht aus den Messpunkten rekonstruieren kann - das wäre aber für
meine Anwendung ein echter Overkill (siehe unten).
Frank: Digitales upsampling würde mein Problem lösen, ich würde aber den
damit verbundenen Rechenaufwand gerne vermeiden. Das VI soll die
Messpunkte in einem Graphen darstellen (kein Problem, sieht dann eben
aus wie oben). Weitere Graphen sollen das Spektrum des Signals anzeigen
(durch FFT realisiert) sowie das 1. und 2. Integral des Signals als
Zeitverlauf oder Spektrum.
Schon beim Spektrum bekomme ich aber sobald die Abtastrate weniger als
4* die gemessene Frequenz ist eine deutlich zu niedrige Amplitude, was
sich meiner Meinung nach aus der Schwankung der Messwerte ergibt. Bei
dem Beispiel von dem das obige Bild stammt (Abtastrate = 2,038*f), ergab
die FFT nur eine Amplitude von 0,68 (Tatsächliche Amplitude, eingespeist
und aus den Maxima des Zeitsignals bestätigt: 1)
Diesen Fehler würde ich gerne durch eine einfache Multiplikation mit
einem Faktor in Abhängigkeit von Abtastrate und gemessener Frequenz
herausrechnen, falls möglich. Da schnelle Messungen mit mehreren
Sensoren (jeweils möglicherweise mit mehreren überlagerten Schwingungen)
geplant sind, wäre der Rechenaufwand für ein digitales Oversampling zu
hoch.
Ok, Du verrätst also nicht, was 1. und 2. Integral ist?
Und Du sagst auch nicht, was Du aus dem Signal berechnen willst?
Es wird schwer weiterzuhelfen.
Denk' mal drüber nach! ;-)
- Frank
Sorry, das 1. Integral ist das Integral x dt (wie von dir gepostet), das
2. Integral ist das 1. Integral nochmal nach dt integriert (Integral
Integral x dt²).
Mit dem Signal will ich persönlich erstmal außer der Darstellung in den
Graphen und den Integralen nichts berechnen - wie in meinem letzten
Beitrag geschrieben ist das Ziel meiner Frage, herauszufinden ob ich die
im Frequenzbereich zu niedrig angezeigte Amplitude (0,68 statt 1 im
Beispiel) anhand Frequenz und Abtastrate herausrechnen kann, statt erst
aufwändig ein digitales Oversampling zu betreiben und dann das daraus
entstehende Signal durch die FFT (fast fourier transformation) zu jagen
;)
Ich denke da gibt es schon mathematisch keine andere Lösung als mehrere
Samples zu verrechnen.
Fsample/2=2500, Fsignal=2450 macht 50Hz zur Nyquistgrenze.
Ist der Quotient Signalfreqenz/Samplefrequenz genau bekannt und konstant
kann man hier vereinfachen:
if (Fsignal<0.25)
phase=.25/Fsignal
else
phase=0.5/(1-Fsignal/0.5);
amplitude=sqrt(sqr(sample[0]) + sqr(sample[phase]))
Fsignal normiert auf Fsample.
Doch bei nicht ganzzahliger "phase" gibt es auch hier Ungenauigkeiten.
Es wäre zumindest ein Interpolationsfilter für die genaue subsample
position des zweiten Summands der Wurzel nötig, dieser bestimmt die
Genauigkeit.
Auch dieses kann vereinfacht werden, durch Annahmen über die Frequenzen.
Bei Fsignal=.49 ist phase=25 ganzzahlig.
Mögliche Fsignal errechnen sich von der ganzzahligen Phase:
Fsignal=-0.5*(0.5/phase-1)
Wie genaus sollts denn werden? Untenstehendes Matlab script macht 60dB
oder so. Das Signal wird gefenstert und zu einem analytischen komplexen
Signal ergänzt (Obere Hälfte der FFT wird 0), dann der abs-Wert aus der
Mitte mal 2.
Das geht auch kontinuierlich mit einem Beobachter für ein leicht
gedämpftes System von genau der Frequenz, da kommt für jeden samplewert
eine Amplitudenschätzung raus.
Cheers
Detlef
clear
fs= 2450;
fa=4994.29;
Ap =2000;
tend=0.12;
tt=linspace(0,tend,fa*tend);
tt=tt(1:512);
sig=Ap*cos(2*pi*fs*tt);
sig=sig.*(hanning(512)).';
sps=fft(sig);
sps(257:end)=0;
anasig=ifft(sps);
plot(abs(anasig))
2*abs(anasig(256))
return
Hi Syndic,
siehst Du, jetzt kommen wir der Sache schon näher.
Also ich würde den RMS berechnen, also so:
Der RMS entspricht dann dem Effektivwert. Bei einem Sinus (!) mußt Du
dann noch mit Wurzel 2 multiplizieren.
Das funktioniert immer, allerdings solltest Du Nyquist einhalten.
Ich würde sagen, so macht man das.
Gruss
- Frank
Frank,
hast du das schonmal experimentell z.B. mit einer TK überprüft? Bei
Signalfrequenz z.B. 0.49 Fs?
Ich hoffe du entwickelst nichts wofür ich einmal Geld ausgebe.
Hallo Dogbert,
also ich habe es gerade mit MATLAB probiert und es funktioniert perfekt.
Es ist aufwandsarm und funktioniert auch bei 0.49 Fs.
Wo siehst Du das Problem, das solltest Du schon erklären und keine
polemischen Sprüche absondern.
Ach, und ja, ich entwickle wahrscheinlich Dinge für Du Geld ausgibst.
(Jedesmal, wenn Du ein Mobilfunknetz benutzt ;-)
Sport frei, Dogbert
- Frank
Also jetzt mal ganz ehrlich Leute,
der Algorithmus mißt die Signalleistung. Bitte nicht schwätzen, sondern
ausprobieren. So und damit es JEDER auch nachvollziehen kann, hier die
paar Zeilen für Matlab:
>> f=0.49;>> t=0:10000;>> x=sin(2*pi*f*t);>> plot(t,x)>> x2=x.*x;>> y=filter([0.01],[1 -0.99],x2);>> plot(sqrt(2*y))
Warum muß man unbedingt zeigen, dass man keine Ahnung hat?
Gruss
- Frank
Hast recht, das war nicht ganz korrekt ausgedrückt. Ich wollte zum
Ausdruck bringen, dass mit den Eingangsdaten vom TO mit nur 1 bis 2
Schwebungsperioden nach dem IFR nicht wirklich das rauskommt, was er
wissen möchte.
Schon klar, wenn das Signal konstant bleibt und die Amplitude
hinreichend langsam schwankt, dann kommt schon was raus.
Aber dein IFR braucht einige 100 Perioden, bis es einen einigermaßen
konstanten Wert erreicht hat. Und dann schwingt das Ding klarerweise
weiter. Was, wenn der TO eigendlich Amplitudenschwankungen von ~ 10 ..
50 Perioden finden möchte?
Mach' mal in der 2. Zeile:
t=0:200;
Liebe Freunde der gepflegten Signalverarbeitung,
ich brauchte ein bißchen Spaß und habe die Amplitude des Signals mit nem
Kalman Filter geschätzt. Das Ergebnis ist im angehängten .jpg zu sehen:
rot das Originalsignal, blau das geschätzte Signal und grün dessen
Amplitude.
Klappt gut.
Vorgehensweise (siehe angehängtes Matlab): Man bastelt sich einen
ungedämpften Schwinger (Systemmatrix A) und läßt ihn schwingen,
Ausgangsgröße ist y, Zustandsgrößen sind x.
Dann tut man so, als sei y der Meßwert und schätzt mit einem Kalman
Filter die Zustandsgrößen xd ('xdach') des Modells (wieder Systemmatrix
A). Dazu benutzt man die Verstärkungen K des Kalmanfilters. Ich nehme
allerdings die stationären ks, mit denen klappt das auch noch. Der
update der K ist nur rauskommentiert, mit dem habe ich die ks berechnet.
In einem 'Produktionscode' reduziert sich das Zeug auf die Berechnung
von:
xds=A*xd;
xd = xds+K*(y(k) -C*xds);
schlichter geht's nimmer.
In dem Zeug ist möglicherweise noch nen Bug: Das negative K ist komisch
und die gschätzten Zustandsgrößen haben auch das falsche Vorzeichen,
trotzdem stimmt natürlich die Amplitude.
Der TO ist ja möglicherweise schon vorm' rauen Ton hier geflüchtet, kann
aber gerne nochmal expliziter werden.
War Spaß jewesen!
math rulez
Cheers
Detlef
clear
%alpha=pi/10;
% Das ist der Drehwinkel des Zeigers
%bißchen weniger als pi pro sample
alpha=(2450/4994.29)*2*pi;
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;
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([abs(x(1,:)+j*(x(2,:)));abs(xxd(1,:)+j*(xxd(2,:)))].')
%plot(interim);
return
Hallo,
also das Kalman-Filter, kennt ja bereits die Frequenz des Signals über
die Variable alpha:
A=Ap*[ cos(alpha) -sin(alpha) ; ...
sin(alpha) cos(alpha) ];
Paßt das denn zum Problem? Kann man davon ausgehen, dass man das weiß?
Wenn dies so ist und das Signal nur einen single tone enthält, dann
könnte man ja auch eine max hold Funktion, die gar keinen Aufwand
bedeutet, verwenden.
Was meint Ihr dazu?
Gruss
- Frank
sinus.png
(6,2 KB)
