ich habe eine Frage, die sich an den Forumbeitrag
http://www.gomatlab.de/pwelch-selber-erstellen-t20937.html
anlehnt, jedoch in eine etwas theoretischere Richtung geht. Es geht mir um den Befehl pwelch im Allgemeinen. Ich hoffe Ihr verzeiht mir, dass ich einen neuen Thread für das Thema als sinnvoll erachte.
Zunächst einmal der ganz leicht abgewandelte Code:
Code:
%% Berechnung des Leistungs-Dichte-Spektrums
Ts=0.001;
fs=1/Ts;
%zwei verschiedene Zeitvektoren t1 und t2 definieren
t1 = 0:Ts:10; % Zeitvektor % Nur gerade Werte verwenden ifmod(length(t1),2)~=0
t1=t1(1:length(t1)-1);
disp('t war ungerade->abgeändert') end
t2 = 0:Ts:5; % Zeitvektor % Nur gerade Werte verwenden ifmod(length(t2),2)~=0
t2=t2(1:length(t2)-1);
disp('t war ungerade->abgeändert') end
iflength(x)==length(Signal_3) subplot(3,1,1); plot(t2,x);
title('Signal');
else subplot(3,1,1); plot(t1,x);
title('Signal');
end
%%
%
N = length(x); % Anzahl Messwerte % Segmenterstellung, Fensterung und Durchführung der FFT
y = x; % Signal
M = 0.1*length(x); % Fensterbreite
win = window(@hamming,M); % Fensterung
OL = 0.5; % Überdeckung der Segmente
numoverlap = M * OL;
k = fix((N-numoverlap)/(M-numoverlap)); % Anzahl der dft Spectren
MY = (M/2)+1; % Anzahl der dft Koeffizienten pro dft Spektrum
Sxx = zeros(MY,1); % der zu füllende Vektor
AL = win' * win; % Normierungskonstante
start=1;
stop=start+M-1;
for i=1:k
YY = (fft(y(start:stop).*win',M))';
Sxx = Sxx + YY(1:MY) .*conj(YY(1:MY))/AL;
start = start + M * (1-OL);
stop = start + M - 1;
end
% weitere Normierung
Sxx = [Sxx(1); 2*Sxx(2:MY-1); Sxx(MY)]; % Gleichanteil und Nyquistfrequenz unverändert
Sxx = Sxx./(fs*k);
Pxx = 10*log10(Sxx+eps); % eps beim logarithmieren nicht vergessen um log10(0) = -Inf zu vermeiden
fx = 0:fs/M:(fs/2); % Frequenzvektor 0...fs/2 [Hz]
Folgende Fragen habe ich:
1. Welche Einheit hat die y-Achse des unteren Plots? Wenn das Signal die Einheit m/s^2 besitzt, müsste das doch (m/s^2)^2 sein oder?
2. Warum gibt der Befehl pwelch nicht genau das aus, was beim händisch errechneten pwelch ausgegeben wird, obwohl die eckigen Klammern eigentlich genau den gewählten Parametern der händischen Version entsprechen (das sind die default Parameter von Matlab)?
3. Wie verhält es sich mit der Fensterung der Signale? Bei einer Änderung der Fensterung, ändern sich die Ergebnisse drastisch. So wird sowohl der Maximalwert des Pxx1 Wertes stark verändert als auch der Verlauf wird beeinflusst. Das macht es natürlich schwierig dieses Ergebnis zur Weiterverarbeitung zu verwenden.
4. Übergreifend zu 3: Generell habe ich festgestellt, dass sich das Ergebnis stark verändert bei Variation aller Parameter. Es stellt sich die Frage, wie die Parameter gewählt werden sollten um plausible Ergebnisse zu bekommen, wobei ich nicht in der Lage bin abzuschätzen was plausibel ist und was nicht. So sollte der Parameter der Abtastung fs vermutlich immer der tatsächlichen Signalabtastung entsprechen (hier fs=1/Ts=1000). Wie verhält es sich mit den anderen Parametern?
Ich bin sehr dankbar für alle Anregungen und Kommentare!
Zu 1) Da eine Leistung gemessen wird und bei der Darstelllung in dB -> dBW/Hz. Ich empfehle die Angabe Magnitude [dB/Hz]...dann kann man nichts falsch machen.
Zu 2) Du setzt den Inputparameter window = []. Das bedeutet, dass der Default verwendet wird. Damit wird das Signal in 8 Segmente geteilt. Beim selbst erechneten PSD wird aber die Fensterbreite so festgelegt...
und somit ergeben sich auch mehr als 8 Segmente. Daher wohl der Unterschied-
Zu 3+4) Gegenfrage...was willst du mit dem PSD eigentlich bestimmen? Wenn man solche Funktionen nutzt, sollte man sich immer vorab mit einem Testbeispiel beschäftigen, bei dem man GENAU weiß was herauskommen soll. Man kann hier keine allgemein gültige Aussage machen, welche Einstellungen zu verwenden sind. Das hängt davon was man ermitteln will und wie das verwendete Signal beschaffen ist (Rauschen etc.)
Danke für die Antwort!
1. Okay ich dachte mir schon, dass da eine Leistung steht. Ob das nun in dB angegeben wird oder nicht spielt ja eigentlich keine Rolle.
2. Das war blöd von mir. Flüchtigkeitsfehler... danke
3 und 4.
Leider habe ich kein Testbeispiel an dem ich das validieren kann. Aber generell geht es um folgendes: Die spektrale Leistungsdichte P eines Signals wird mit einer Funktion f die im Frequenzbereich vorliegt mulitpliziert und das Produkt wird im Frequenzbereich integriert.
Es ergibt sich also am Ende ein skalar M, welches zur Auswertung herangezogen wird.
Das eigentliche Signal (aus dem die Leistungsdichte ermittelt wird), schwingt dabei nach plötzlicher Anregung mit einer Amplitude von ca. plus minus 3 und einer Frequenz von bis zu 10Hz. Nach zwei bis drei Perioden ist die Schwingung wieder abgeklungen. Generell liegt permanent ein Grundrauschen mit einer Amplitude von plus minus 0.3 vor.
Obwohl ich mich schon etwas in das Thema eigelesen habe, habe ich nun Schwierigkeiten in der Funktion pwelch die richtigen Parameter zu setzen. Sobald die Fensterung oder andere Parameter verändert werden, ändert sich das P sehr stark... da bräuchte ich noch ein paar Tipps
Ich kann mich da leider nur wiederholen. Allgemein gültige Aussagen zu den richtigen Einstellungen gibt es nicht. Es hängt vom Signal ab und falls ein System untersucht werden soll, ob z.B. ein MA oder AR Prozess vorliegt.
Bei der Schätzung der Leistungsdichte nach Welch wird aber immer ein Hamming Fenster verwendet. Bei einem Rechteck-Fenster (sozusagen der Sonderfall zu Welch) heißt es dann z.B. Bartlett-Methode.
Wählt man die zu mittelnden Periodogramme (Anzahl Fenster) zu hoch, ist die Schätzung gegenüber der tatsächlichen Leistungsdichte zu gering. Außerdem wird dadurch die Frequenzauflösung schlechter. Sind die Fenster dagegen zu lang gewählt, liegt die Schätzung über dem tatsächlichen Wert.
Du wirst bei dem Thema nicht drum herumkommen Lektüre zu lesen. Z.B. Digitale Signalverarbeitung von K.D. Kammeyer mit Matlab Übungen.
Im Übrigen lagst du mit deiner Annahme der Einheit doch richtig. Wird z.B. eine Schwinungungsanlayse eines mechanischen Systems durchgeführt, spricht man von einer spektralen Beschleunigungsdichte.
Okay, ich werde dann versuchen mich noch ein wenig mehr einzulesen.
Danke für die Hilfe!!
Einstellungen und Berechtigungen
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