Verfasst am: 29.06.2011, 12:34
Titel: Grundwelle eines Signals berechnen
Hallo Forenmitglieder!
Ich bin neu im Gebiet der Signalverarbeitung und insbesondere in der Anwendung von MATLAB/Simulink.
Mich beschäftigt nun folgendes Problem:
Ich habe ein sinusförmiges Signal (aber verrauscht) welches aus N diskreten Zeit- t(0...T) und den zugehörigen Funktionswerten x(t(0...T)) besteht.
Gemessen wurde also genau eine Periode T des Signals.
Sowohl die Zeit- als auch die Funktionswerte befinden sich als Vektor in MATLAB-Variablen.
Ich würde nun gerne die Grundwelle dieses Signals extrahieren um das Signal damit in erster Näherung approximieren zu können.
Weiss jemand, wie ich dies in MATLAB realisieren kann?
1. das signal mit einer gewissen Eckfrequenz Tiefpass-Filtern
2. eine FFT rechnen bzw. das Spektrum analysieren
3. eine Funktion z.B. Sinus in das signal fitten
was hätten wir denn gerne?
Problem bei
1. wahrscheinlich bekommt man einen störenden Phasenfehler im Signal
2. die Qualität hängt ab von der Anzahl an Messpunkten
deshalb:
kennt man das Signal so ungefähr, so kann man ein curve fitting durchführen und das Signal ideal "erzeugen"
Hallo!
Erstmal vielen Dank für eure interessanten Antworten.
@Idefix_1024:
Die Curve Fitting Toolbox besitze ich leider nicht. Möglichkeit 2 klingt jedoch sehr interessant. Messpunkte sollten genügend vorhanden sein.
Kannst du vielleicht etwas mehr dazu sagen bzw. sogar Beispielcode geben?
Die FFT zerlegt dein Signal in die vorhandenen Einzelsignale. Auf der X-Achse hast du die Frequenz und auf der Y-Achse die Amplitude. Damit wirst du auch deine Grundschwingung des Signals (als Peak) sehen und kannst dann die dazugehörige Freq. ablesen.
hier ein Beispiel bei dem klar werden sollte wie das gemeint war
Code:
%% Signal
dt = 1/8000;
zeit = 0:dt:0.1;
f1 = 100;
f2 = 500;
Ampl1 = 10;
Ampl2 = 5;
offset = 2;
signal = offset+(Ampl1*cos(2*pi*f1*zeit)+Ampl2*cos(2*pi*f2*zeit));
Ts = diff(zeit(1:2)); % Abtastzeit
N = length(zeit); % Länge des Datenvektors
f = [0:floor((N-1)/2)] / (N*Ts); % Frequenzvektor für Spektrum-Plot % Fouriertransformierte
SIGNAL = fft(signal); % Fouriertransformation
SIGNAL = SIGNAL / N; % Normierung
SIGNAL = [SIGNAL(1)2*SIGNAL(2:floor((N-1) / 2) + 1)]; % begrenzen auf f_max
ergebnis = [SIGNAL; f];
%% Plot der Berechnungen figure() subplot(211) plot(zeit,signal,'r','LineWidth',2);
xlabel('Zeit [s]');
ylabel('Signal');
grid on
subplot(212) plot(f,abs(SIGNAL),'rx','LineWidth',2,'Markersize',8);
xlabel('Frequenz [Hz]');
ylabel('Absolutwert');
grid on
Super!
Vielen Dank für den Code. Anhand der Plots verstehe ich nun was gemeint ist.
Drei Verständnisfragen jedoch zu deinem Skript:
1.) Wie kommst du auf den Frequenzvektor für den Spektrumplot?
2.) Warum machst du das mit der Normierung
3.) Warum ist die Begrenzung auf f_max erforderlich und wie kommst du darauf
Das Signal besteht doch hauptsächlich aus einer Grundschwingung und irgendwie geartetem Rauschen. da T bekannt ist, sollte auch f kein Problem sein, ohne den Zeitbereich verlassen zu müssen... ich seh das eher als Scherzfrage eines Profs an, der sehen will, wer alles mit FFT loslegt
bei dem Satz
Idefix_1024 hat Folgendes geschrieben:
ein Sinus-Signal besteht aber nunmal nicht nur aus der Frequenz...
überleg ich, ob ich das Problem richtig gesehen habe... warhscheinlich nicht...
Wie oben beschrieben, stell ich mir das Signal mit einer großen Grundschwingung vor, und dessen Frequenz meinte ich. Dass ist die, die bestimmt werden soll. und von der wurde eine Periode aufgenommen.
Wo ist jetzt mein Denkfeher?
_________________
Ich hasse es wenn die Leute Fragen stellen, man dann versucht sich Mühe zu geben, und diejenigen ihren Thread nie wieder besuchen...
Du hast da prinzipiell keinen Denkfehler. Ist wirklich genau eine Periode der Grundschwingung gemessen, kann man die Frequenz leicht bestimmen. Bei der Amplitude und Phase sieht es da schon schwieriger aus...das kommt halt auf die überlagerte Störungen an. Auf dem Bild ist das schon schwerer zu erkennen..ich habe hier folgende Daten aus dem Skript von Idefix verwendet (da war das Signal nämlich mehr als eine Periode). In diesem Bsp. ist f2 ja noch günstig gewählt, da es ein Vielfaches von f1 ist. Ist das nicht der Fall, wird's noch schwieriger.
Code:
% Signal
dt = 1/8000;
zeit = 0:dt:0.01; % genau 1 Periode der Grundschwingung f1
f1 = 100;
f2 = 500;
Ampl1 = 1;
Ampl2 = 5;
offset = 2;
Also so wie ich das sehe besteht ein sinusförmiges Signal IMMER aus einer Amplitude, einer Phase und einer Frequenz.
Möchte man also die Grundwelle ideal darstellen, so benötigt man alle drei Parameter. Da genügt es eher selten nur die Frequenz zu kennen denke ich.
Sowohl Amplitude, als auch Phase lassen sich mit einer FFT gewinnen (vorausgesetzt man kann mit einer gewissen Unschärfe, die von der Anzahl der Messpunkte abhängt leben).
Es geht aus meiner Sicht nicht nur um die Frequenz. Die könnte man aus jedem Zeitverlauf durch "Hinschauen" bestimmen. Da hast du Recht Andy386!
Zu den drei Fragen kann ich leider spontan nur sagen, dass sich auf Grund des Abtasttheorems (pseudowissenschaftl. Kurzfassung: ein Signal mit einer Frequenz f muss man mit einer doppelt so hohen Frequenz abtasten um es noch richtig darstellen zu können) der maximal zu analysierende Frequenzbereich genau zur Hälfte ergibt. Im Beispiel sind die Messdaten mit 8 kHz abgetastet worden. Deshalb macht nur eine FFT bis 4 kHz Sinn.
Die Normierung ergibt sich aus der mathematischen Berechnung der FFT und ist Teil des Algorithmus. Der Frequenzvektor eigentlich ebenso. Die FFT liefert zu jeder Stützstelle im Spektrum eine komplexe Zahl. In dieser komplexen Zahl steckt Phase und Betrag zu dieser Frequenz.
Vielleicht hilft Dir auch Google zum weiteren Verständnis weiter. Wie genau die FFT funktioniert ist kaum auf zwei Zeilen zu erklären ;-)
wenn Du denkst, dass eine Filterung besser ist, dann schau Dir einfach den direkten Vergleich an...
Code:
%% Signal
dt = 1/8000;
zeit = 0:dt:0.1;
f1 = 100;
f2 = 500;
Ampl1 = 10;
Ampl2 = 5;
offset = 2;
signal = offset+(Ampl1*cos(2*pi*f1*zeit)+Ampl2*cos(2*pi*f2*zeit));
Ts = diff(zeit(1:2)); % Abtastzeit
N = length(zeit); % Länge des Datenvektors
f = [0:floor((N-1)/2)] / (N*Ts); % Frequenzvektor für Spektrum-Plot % Fouriertransformierte
SIGNAL = fft(signal); % Fouriertransformation
SIGNAL = SIGNAL / N; % Normierung
SIGNAL = [SIGNAL(1)2*SIGNAL(2:floor((N-1) / 2) + 1)]; % begrenzen auf f_max
% Fourierkoeffizienten a und b bestimmen for i=1:1:length(SIGNAL)
a(i) = 2*real(SIGNAL(i))/N;
b(i) = -2*imag(SIGNAL(i))/N;
end;
a(1) = a(1)/2;
% Frequenzvektor sortieren so dass später das Signal aus den Hauptfrequenzanteilen % wieder erzeugt bzw. approximiert werden kann [wert position] = sort(abs(SIGNAL));
for i=1:1:length(f)
f_(i) = f(position(length(position)-i+1));
a_(i) = a(position(length(position)-i+1)) .*N;
b_(i) = b(position(length(position)-i+1)) .*N;
end
%% Signalsynthese auf Grund der FFT
signal_FFT = Fourier_Signal_Syth(signal, zeit, 1);
%% Plot der Berechnungen figure(1) subplot(211) plot(zeit,signal,'r','LineWidth',2);
xlabel('Zeit [s]');
ylabel('Signal');
grid on
subplot(212) plot(f,abs(SIGNAL),'rx','LineWidth',2,'Markersize',8);
xlabel('Frequenz [Hz]');
ylabel('Absolutwert');
grid on
figure(2) hold on
plot(zeit,signal,'r','LineWidth',2);
plot(zeit,signal_FFT,'b','LineWidth',2);
plot(zeit,signal_filt,'m','LineWidth',2);
hold off
xlabel('Zeit [s]');
ylabel('Signal');
grid on
und was fällt auf?
Mit diesem Skript kann man den Unterschied sofort sehen, den ich mit Phasenfehler gemeint habe. Außerdem wird das Problem des Initialisieren des Filters deutlich. Ein Filter muss immer erst einschwingen bevor es brauchbare Ergebnisse liefert.
Ja...hast Recht. Das Filter ist hier ungünstig, wenn wirklich nur eine Periode vorhanden ist. In diesem Fall war das kein guter Tipp, da sich das MW Filter nicht rechtzeitig einschwingt. Aber bei einer längeren Messzeit ist ein Filter sicherlich genauso hilfreich und wesentlich einfacher zu implementieren als eine FFT...vor allem dann wenn man sich damit nicht auskennt
Du kannst Beiträge in dieses Forum schreiben. Du kannst auf Beiträge in diesem Forum antworten. Du kannst deine Beiträge in diesem Forum nicht bearbeiten. Du kannst deine Beiträge in diesem Forum nicht löschen. Du kannst an Umfragen in diesem Forum nicht mitmachen. Du kannst Dateien in diesem Forum posten Du kannst Dateien in diesem Forum herunterladen
MATLAB, Simulink, Stateflow, Handle Graphics, Real-Time Workshop, SimBiology, SimHydraulics, SimEvents, and xPC TargetBox are registered trademarks and The MathWorks, the L-shaped membrane logo, and Embedded MATLAB are trademarks of The MathWorks, Inc.
RSS