mikrocontroller.net

Forum: Digitale Signalverarbeitung / DSP Sinus mit 2f < Abtastrate < 4*f: Schwankung herausrechnen


Autor: Syndic (Gast)
Datum:
Angehängte Dateien:

Bewertung
0 lesenswert
nicht lesenswert
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.

Autor: Frank (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
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

Autor: Mark Brandis (markbrandis)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
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...

Wenn dem nicht so ist, dann kann Deine Aussage mit Abtastrate > 
2*Signalfrequenz wohl nicht stimmen.

Autor: Mark Brandis (markbrandis)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
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.9Cb...

Autor: Syndic (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
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.

Autor: Frank (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
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

Autor: Syndic (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
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 
;)

Autor: Dogbert (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
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)

Autor: Detlef _a (detlef_a)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
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

Autor: Zwie Blum (zwieblum)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
@Syndic:

Das schaut aber sehr verdächtig nach Schwebung aus. Was sagt denn dein 
Oszi zum Eingangssignal?

Autor: Frank (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
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

Autor: Dogbert (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
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.

Autor: Frank (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
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

Autor: Zwie Blum (zwieblum)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Dein IFR wirft dir die Hüllkurve der Schwebung aus. Das ist aber nicht 
das, was der TO berechnen will.

Autor: Helmut Lenzen (helmi1)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Dir ist bekannt das die Amplitude nach einer sinc(x) = sin(x)/x Funktion 
geht?

http://de.wikipedia.org/wiki/Nyquist-Shannon-Abtasttheorem

Autor: Frank (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
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

Autor: Zwie Blum (zwieblum)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Deine Annahme der Eingangsdaten ist falsch.

Autor: Frank (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Nein, die ist genau so, wie oben angegeben.

Hey, Zwie Blum, Du mußt schon ein wenig konkreter werden und nicht nur 
schlaumeinern.

- Frank

Autor: Zwie Blum (zwieblum)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
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;

Autor: Detlef _a (detlef_a)
Datum:
Angehängte Dateien:

Bewertung
0 lesenswert
nicht lesenswert
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

Autor: Zwie Blum (zwieblum)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Fesch :-)

Autor: Frank (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
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

Autor: Dogbert (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Detlefs Kalman Filter muss doch sogar die Phase des Signals kennen, denn 
die ist immer 0 wenn ich nicht irre.

Antwort schreiben

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

Formatierung (mehr Informationen...)

  • [c]C-Code[/c]
  • [avrasm]AVR-Assembler-Code[/avrasm]
  • [code]Code in anderen Sprachen, ASCII-Zeichnungen[/code]
  • [math]Formel in LaTeX-Syntax[/math]
  • [[Titel]] - Link zu Artikel
  • Verweis auf anderen Beitrag einfügen: Rechtsklick auf Beitragstitel,
    "Adresse kopieren", und in den Text einfügen




Bild automatisch verkleinern, falls nötig
Bitte das JPG-Format nur für Fotos und Scans verwenden!
Zeichnungen und Screenshots im PNG- oder
GIF-Format hochladen. Siehe Bildformate.
Hinweis: der ursprüngliche Beitrag ist mehr als 6 Monate alt.
Bitte hier nur auf die ursprüngliche Frage antworten,
für neue Fragen einen neuen Beitrag erstellen.

Mit dem Abschicken bestätigst du, die Nutzungsbedingungen anzuerkennen.