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.
Angehängte Dateien:
-
sinus.png
6,2 KB
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:
Beschreib nochmal genau! Gruss - Frank
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
@Syndic: Das schaut aber sehr verdächtig nach Schwebung aus. Was sagt denn dein Oszi zum Eingangssignal?
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
Dein IFR wirft dir die Hüllkurve der Schwebung aus. Das ist aber nicht das, was der TO berechnen will.
Dir ist bekannt das die Amplitude nach einer sinc(x) = sin(x)/x Funktion geht? http://de.wikipedia.org/wiki/Nyquist-Shannon-Abtasttheorem
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
Deine Annahme der Eingangsdaten ist falsch.
Nein, die ist genau so, wie oben angegeben. Hey, Zwie Blum, Du mußt schon ein wenig konkreter werden und nicht nur schlaumeinern. - 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;
Angehängte Dateien:
-
kal.jpg
67 KB
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
Detlefs Kalman Filter muss doch sogar die Phase des Signals kennen, denn die ist immer 0 wenn ich nicht irre.
Thread beobachten
Seitenaufteilung abschalten
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.