kennt sich jemand mit dem Algorithmus von pwelch aus? Ich versuche gerade den nachzuahmen allerdings scheitere ich kläglich. Mir fehlt die Signaprocessing Toolbox und somit kann ich den nicht nutzen.
Was ich bisher herausgefunden habe ist, dass man im Grunde die Zeitsignale in Segmente aufteilt, darüber eine FFT laufen lässt, diese quadriert und normiert und Anschließend alle Segmente wieder aufsummiert. Das ist mein derzeitiges Wisen, allerdings bin ich mir nicht sicher ob das wirklich so geht!
Du hast noch die Fensterung jedes Segments vor der FFT vergessen. Bei pwelch wird standardmäßig ein Hamming window verwendet.
Poste doch bitte mal deinen Code, damit man nicht von Null anfangen muss. Bis zur Skalierung der transform. Segmente sollte ja alles zu dem Code des Spektrogramms identisch sein, den du ja schon hast
%% Berechnung des Leistungs-Dichte-Spektrums (nach Welch) mit einem Sinus
f = 1000; % Frequenz des zu erzeugenden Sinus
t = 0:1/16383:1; % Zeitvektor
x = sin(2*pi*f*t); % Sinusvektor
N = length(x); % Anzahl Messwerte
% Segmenterstellung, Fensterung und Durchführung der FFT
fs = length(t);
y = x; % Signal
M = fs/16; % Fensterbreite
win = Fensterung(M,'hamming'); % 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 + abs(YY(1:MY)).^2/AL; % aufsummieren,quadrieren und normieren der einzelnen Periodogramme % weitere Möglichkeit
%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; % siehe Matlab-Doku zu pwelch
Sxx = Sxx./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]
vielen, vielen, vielen Dank. Das ist ja wunderbar.
Jetzt muss ich das noch ein wenig umschreiben, damit ich zu jedem Segment den Schwerpunkt bekomme. Dazu muss ich ja im Grunde den Sxx erweitern, so dass alle FFT´serstmal komplett gespeichert werden ( in einer Matrix ) und da alle Umrechnungen dran vorgenommen werden. Die Aufsummierung müsste ich ja auch zum Schluss erledigen können!
Sollte mal wieder ein Denkfehler vorhanden sein, korrigier mich doch bitte kurz.
solltest du öfter Algorithmen aus der Signal Processing Toolbox brauchen, würde ich mir die Anschaffung dieser Toolbox überlegen. Warum viel Zeit und Energie investieren, um etwas zu schaffen, dass es schon gibt?
Eine Sache hatte ich noch vergessen zu erwähnen. In deinem Code fehlt die Abfrage, ob die Anzahl N von Messwerten gerade oder ungerade ist. Beim Spektrogramm wurde ja immer von gerader Anzahl ausgegangen und wenn nötig sonst der letzte Messwert abgeschnitten. Das solltest du hier auch einbauen. Falls du immer alle Messwerte verwenden möchtest, musst du die Fälle z.B. bei der Normierung getrennt abarbeiten.
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