clear all; close all; clc %% Einlesen der einzelnen 10 Messdaten % ------------------------------------ % Daten mit dem Iindex M stehen für den defekten Lüfter (er ist % absichtlich manipuliert worden, deswegen M für manipuliert)! % Die Kanalbezeichnungen stammen aus den vier Eingängen des PicoScope % Oszilloskop (Kanal A, Kanal B, Kanal C und Kanal D) Datei1 = load('D:\Desktop\Messdaten\I_PWM100_Speed5720_hochlauf.mat\I_PWM100_Speed5720_hochlauf_01.mat'); Datei1c = struct2cell(Datei1); KanalA_1 = Datei1c{6,1}; KanalB_1 = Datei1c{7,1}; KanalC_1 = Datei1c{8,1}; KanalD_1 = Datei1c{9,1}; Datei1M = load('D:\Desktop\Messdaten\I_M_PWM100_Speed5830_hochlauf.mat\I_M_PWM100_Speed5830_hochlauf_01.mat'); Datei1cM = struct2cell(Datei1M); KanalA_1M = Datei1cM{6,1}; KanalB_1M = Datei1cM{7,1}; KanalC_1M = Datei1cM{8,1}; KanalD_1M = Datei1cM{9,1}; %% Messdaten in eine Datei zusammenfügen % -------------------------------------- KanalA_kpl = vertcat(KanalA_1); KanalAM_kpl = vertcat(KanalA_1M); KanalB_kpl = vertcat(KanalB_1); KanalBM_kpl = vertcat(KanalB_1M); KanalC_kpl = vertcat(KanalC_1); KanalCM_kpl = vertcat(KanalC_1M); KanalD_kpl = vertcat(KanalD_1); KanalDM_kpl = vertcat(KanalD_1M); %% Variablenzuweisung % ------------------- % Rohdaten in mV % -------------- x = KanalA_kpl;% Daten in x-Achse - standard Lüfter xM = KanalAM_kpl;% Daten in x-Achse - defekter Lüfter y = KanalB_kpl;% Daten in y-Achse - standard Lüfter yM = KanalBM_kpl;% Daten in y-Achse - defekter Lüfter z = KanalC_kpl;% Daten in z-Achse - standard Lüfter zM = KanalCM_kpl;% Daten in z-Achse - defekter Lüfter mic = KanalD_kpl;% Mikrofon - standard Lüfter micM = KanalDM_kpl;% Mikrofon - defekter Lüfter % Beschleunnigungssensor % ---------------------- % Daten in m/s^2 % -------------- xa = x/1.011;% Umrechnungsfaktor lt. Kalibrierzertifikat (1,011 mV/(m/s^2)) xaM = xM/1.011;% Umrechnungsfaktor lt. Kalibrierzertifikat (1,011 mV/(m/s^2)) ya = y/1.006;% Umrechnungsfaktor lt. Kalibrierzertifikat (1,006 mV/(m/s^2)) yaM = yM/1.006;% Umrechnungsfaktor lt. Kalibrierzertifikat (1,006 mV/(m/s^2)) za = z/0.988;% Umrechnungsfaktor lt. Kalibrierzertifikat (0,988 mV/(m/s^2)) zaM = zM/0.988;% Umrechnungsfaktor lt. Kalibrierzertifikat (0,988 mV/(m/s^2)) % Mikrofon % -------- % Daten in Pa % ----------- micp = mic/50.2;% Umrechnungsfaktor lt. Kalibrierzertifikat (50,2 mV/Pa) micpM = micM/50.2;% Umrechnungsfaktor lt. Kalibrierzertifikat (50,2 mV/Pa) % Parameter % --------- fs = 50000;% Anzahl Messpunkte pro Sekunde (Abtastfrequenz) n = length(x);% Anzahl der gesamten Messpunkte x-Achse nfft = 16384;% nfft =2^14, Sampels für die Frequenzauflösung fn = fs/2;% Nyquist-Frequenz t = (0:n-1)/fs;% Messzeit dt = 1/fs;% Zeitschritt df = fs/nfft;% Frequenzauflösung x_fn = 0 : df : fn-df;% Nyquist-Frequenzvektor x_fs = 0 : df : fs-df;% Abtast-Frequenzvektor bw_a = 1e-6;% Bezugswert für Beschleunigung nach DIN EN ISO 1683 in mikro_m/s^2 bw_p = 1e-6;% Bezugswert für Schalldruck nach DIN EN ISO 1683 in mikro_Pa %% Beschleunigungssensor - Darstellung im Zeitbereich % --------------------------------------------------- % Achsenskalierung (wird für die graphische Darstellung benötigt) % ---------------- % Ermittlung der kleinsten Amplitude if (min(xa) <= min(ya)) min_y = min(xa); elseif (min(ya) <= min(za)) min_y = min(ya); elseif (min(za) <= min(xaM)) min_y = min(za); elseif (min(xaM) <= min(yaM)) min_y = min(xaM); elseif (min(yaM) <= min(zaM)) min_y = min(yaM); else min_y = min(zaM); end % Ermittlung der größten Amplitude if (max(xa) >= max(ya)) max_y = max(xa); elseif (max(ya) >= max(za)) max_y = max(ya); elseif (max(za) >= max(xaM)) max_y = max(za); elseif (max(xaM) >= max(yaM)) max_y = max(xaM); elseif (max(yaM) >= max(zaM)) max_y = max(yaM); else max_y = max(zaM); end % Graphische Darstellung % ---------------------- figure('Name','Zeitbereich - Körperschall - standard Lüfter vs. defekter Lüfter - PWM 100 %', ... 'NumberTitle','Off');% Beschriftung Hauptüberschrift plot(t, xa, 'k','LineWidth',1.0);% plot Daten x-Achse - standard Lüfter hold on;% in gleiches Bild plotten plot(t, xaM, 'r','LineWidth',1.0);% plot Daten x-Achse - defekter Lüfter plot(t, ya, 'k','LineWidth',1.0);% plot Daten y-Achse - standard Lüfter plot(t, yaM, 'r','LineWidth',1.0);% plot Daten y-Achse - defekter Lüfter plot(t, za, 'k','LineWidth',1.0);% plot Daten z-Achse - standard Lüfter plot(t, zaM, 'r','LineWidth',1.0);% plot Daten z-Achse - defekter Lüfter set(gca,'FontSize',12);% Schriftgröße set(gcf,'color','w');% weißer Hintergrund axis([0 max(t) min_y max_y]);% Skalierung der x- und y-Achse title('Zeitbereich');% Titel xlabel('Zeit in s');% Beschriftung x-Achse ylabel('Amplitude in m/s^2');% Beschriftung y-Achse legend('standard Lüfter', 'defekter Lüfter','Location','best','LineWidth',1.0);% Legende grid on;% Gitternetzlinie hold off;% nicht mehr in gleiches Bild plotten %% Mikrofon - Darstellung im Zeitbereich % -------------------------------------- % Achsenskalierung (wird für die graphische Darstellung benötigt) % ---------------- % Ermittlung der kleinsten Amplitude if (min(micp) <= min(micpM)) min_yp = min(micp); else min_yp = min(micpM); end % Ermittlung der größten Amplitude if (max(micp) >= max(micpM)) max_yp = max(micp); else max_yp = max(micpM); end % Graphische Darstellung % ---------------------- figure('Name','Zeitbereich - Luftschall - standard Lüfter vs. defekter Lüfter - PWM 100 %', ... 'NumberTitle','Off');% Beschriftung Hauptüberschrift plot(t, micp, 'k','LineWidth',1.0);% plot Daten Mikrofon - standard Lüfter hold on; plot(t, micpM, 'r','LineWidth',1.0);% plot Daten Mikrofon - defekter Lüfter set(gca,'FontSize',12);% Schriftgröße set(gcf,'color','w');% weißer Hintergrund axis([0 max(t) min_yp max_yp]);% Skalierung der x- und y-Achse title('Zeitbereich');% Titel xlabel('Zeit in s');% Beschriftung x-Achse ylabel('Amplitude in Pa');% Beschriftung y-Achse legend('standard Lüfter', 'defekter Lüfter','Location','best','LineWidth',1.0);% Legende, Ausrichtung grid on;% Gitternetzlinie hold off;% nicht mehr in gleiches Bild plotten %% Fensterfunktion % ---------------- % x- Achse - standard Lüfter % -------------------------- k_x = floor(n/nfft);% k Segmente der Länge nfft transformieren H_x = zeros(nfft,1);% Zeropadding, fehlende Daten werden mit Nullen aufgefüllt Hpos0_x = zeros(nfft,k_x); m_x = 1; win_x = hann(nfft); for i_x=1:k_x % Fensterung mit Hann x_win = xa(m_x:nfft+m_x-1).*win_x; H_x = fft(x_win, nfft); % Betrag bilden und positiven Frequenzbereich speichern H_pos0_x(1:(nfft/2)+1,i_x) = abs(H_x(1:(nfft/2)+1)); % nächstes Segment m_x = m_x + nfft; end % x-Achse - defekter Lüfter % ------------------------- k_xM = floor(n/nfft);% k Segmente der Länge nfft transformieren H_xM = zeros(nfft,1);% Zeropadding, fehlende Daten werden mit Nullen aufgefüllt Hpos0_xM = zeros(nfft,k_xM); m_xM = 1; win_xM = hann(nfft); for i_xM=1:k_xM % Fensterung mit Hann xM_win = xaM(m_xM:nfft+m_xM-1).*win_xM; H_xM = fft(xM_win, nfft); % Betrag bilden und positiven Frequenzbereich speichern H_pos0_xM(1:(nfft/2)+1,i_xM) = abs(H_xM(1:(nfft/2)+1)); % nächstes Segment m_xM = m_xM + nfft; end % y-Achse - standard Lüfter % ------------------------- k_y = floor(n/nfft);% k Segmente der Länge nfft transformieren H_y = zeros(nfft,1);% Zeropadding, fehlende Daten werden mit Nullen aufgefüllt Hpos0_y = zeros(nfft,k_y); m_y = 1; win_y = hann(nfft); for i_y=1:k_y % Fensterung mit Hann y_win = ya(m_y:nfft+m_y-1).*win_y;% Multiplikation von Signal und Fenster H_y = fft(y_win, nfft); % Betrag bilden und positiven Frequenzbereich speichern H_pos0_ya(1:(nfft/2)+1,i_y) = abs(H_y(1:(nfft/2)+1)); % nächstes Segment m_y = m_y + nfft; end % y-Achse - defekter Lüfter % ------------------------- k_yM = floor(n/nfft);% k Segmente der Länge nfft transformieren H_yM = zeros(nfft,1);% Zeropadding, fehlende Daten werden mit Nullen aufgefüllt Hpos0_yM = zeros(nfft,k_yM); m_yM = 1; win_yM = hann(nfft); for i_yM=1:k_yM % Fensterung mit Hann yM_win = yaM(m_yM:nfft+m_yM-1).*win_yM;% Multiplikation von Signal und Fenster H_yM = fft(yM_win, nfft); % Betrag bilden und positiven Frequenzbereich speichern H_pos0_yaM(1:(nfft/2)+1,i_yM) = abs(H_yM(1:(nfft/2)+1)); % nächstes Segment m_yM = m_yM + nfft; end % z-Achse - standard Lüfter % ------------------------- k_z = floor(n/nfft);% k Segmente der Länge nfft transformieren H_z = zeros(nfft,1);% Zeropadding, fehlende Daten werden mit Nullen aufgefüllt Hpos0_z = zeros(nfft,k_z); m_z = 1; win_z = hann(nfft); for i_z=1:k_z % Fensterung mit Hann z_win = za(m_z:nfft+m_z-1).*win_z;% Multiplikation von Signal und Fenster H_z = fft(z_win, nfft); % Betrag bilden und positiven Frequenzbereich speichern H_pos0_za(1:(nfft/2)+1,i_z) = abs(H_z(1:(nfft/2)+1)); % nächstes Segment m_z = m_z + nfft; end % z-Achse - defekter Lüfter % ------------------------- k_zM = floor(n/nfft);% k Segmente der Länge nfft transformieren H_zM = zeros(nfft,1);% Zeropadding, fehlende Daten werden mit Nullen aufgefüllt Hpos0_zM = zeros(nfft,k_zM); m_zM = 1; win_zM = hann(nfft); for i_zM=1:k_zM % Fensterung mit Hann zM_win = zaM(m_zM:nfft+m_zM-1).*win_zM;% Multiplikation von Signal und Fenster H_zM = fft(zM_win, nfft); % Betrag bilden und positiven Frequenzbereich speichern H_pos0_zaM(1:(nfft/2)+1,i_zM) = abs(H_zM(1:(nfft/2)+1)); % nächstes Segment m_zM = m_zM + nfft; end % Mikrofon - standard Lüfter % -------------------------- k_mic = floor(n/nfft);% k Segmente der Länge nfft transformieren H_mic = zeros(nfft,1);% Zeropadding, fehlende Daten werden mit Nullen aufgefüllt Hpos0_mic = zeros(nfft,k_mic); m_mic = 1; win_mic = hann(nfft); for i_mic=1:k_mic % Fensterung mit Hann mic_win = micp(m_mic:nfft+m_mic-1).*win_mic;% Multiplikation von Signal und Fenster H_z = fft(mic_win, nfft); % Betrag bilden und positiven Frequenzbereich speichern H_pos0_micp(1:(nfft/2)+1,i_mic) = abs(H_mic(1:(nfft/2)+1)); % nächstes Segment m_mic = m_mic + nfft; end % Mikrofon - defekter Lüfter % -------------------------- k_micM = floor(n/nfft);% k Segmente der Länge nfft transformieren H_micM = zeros(nfft,1);% Zeropadding, fehlende Daten werden mit Nullen aufgefüllt Hpos0_micM = zeros(nfft,k_micM); m_micM = 1; win_micM = hann(nfft); for i_micM=1:k_micM % Fensterung mit Hann micM_win = micpM(m_micM:nfft+m_micM-1).*win_micM;% Multiplikation von Signal und Fenster H_micM = fft(micM_win, nfft); % Betrag bilden und positiven Frequenzbereich speichern H_pos0_micpM(1:(nfft/2)+1,i_micM) = abs(H_micM(1:(nfft/2)+1)); % nächstes Segment m_micM = m_micM + nfft; end %% Beschleunigungssensor - Darstellung im Bildbereich % --------------------------------------------------- % x-Achse - standard Lüfter % ------------------------- X = fft(x_win, nfft);% FFT mit Fensterfunktion Xabs = abs(X);% Absolutzahlen Xfin = fftshift(Xabs/nfft);% Verschiebung der Null-Frequenz-Komponente % zur Mitte des Spektrums % x-Achse - defekter Lüfter % ------------------------- XM = fft(xM_win, nfft);% FFT mit Fensterfunktion XMabs = abs(XM);% Absolutzahlen XMfin = fftshift(XMabs/nfft);% Verschiebung der Null-Frequenz-Komponente % zur Mitte des Spektrums % y-Achse - standard Lüfter % ------------------------- Y = fft(y_win, nfft);% FFT mit Fensterfunktion Yabs = abs(Y);% Absolutzahlen Yfin = fftshift(Yabs/nfft);% Verschiebung der Null-Frequenz-Komponente % zur Mitte des Spektrums % y-Achse - defekter Lüfter % ------------------------- YM = fft(yM_win, nfft);% FFT mit Fensterfunktion YMabs = abs(YM);% Absolutzahlen YMfin = fftshift(YMabs/nfft);% Verschiebung der Null-Frequenz-Komponente % zur Mitte des Spektrums % z-Achse - standard Lüfter % ------------------------- Z = fft(z_win, nfft);% FFT mit Fensterfunktion Zabs = abs(Z);% Absolutzahlen Zfin = fftshift(Zabs/nfft);% Verschiebung der Null-Frequenz-Komponente % zur Mitte des Spektrums % z-Achse - defekter Lüfter % ------------------------- ZM = fft(zM_win, nfft);% FFT mit Fensterfunktion ZMabs = abs(ZM);% Absolutzahlen ZMfin = fftshift(ZMabs/nfft);% Verschiebung der Null-Frequenz-Komponente % zur Mitte des Spektrums % Graphische Darstellungen % ------------------------ % Amplitude in m/s^2 % ------------------ figure('Name','Bildbereich - Körperschall - standard Lüfter vs. defekter Lüfter - PWM 100 %', ... 'NumberTitle','Off');% Beschriftung Hauptüberschrift stem(x_fs-fn, Xfin, 'k.-','LineWidth',3.0);% hold on; stem(x_fs-fn, XMfin, 'r.-','LineWidth',1.0);% stem(x_fs-fn, Yfin, 'k.-','LineWidth',3.0);% stem(x_fs-fn, YMfin, 'r.-','LineWidth',1.0);% stem(x_fs-fn, Zfin, 'k.-','LineWidth',3.0);% stem(x_fs-fn, ZMfin, 'r.-','LineWidth',1.0);% set(gca,'FontSize',12,'XScale','log');% Schriftgröße, logarithmische x-Achse, Skalierung x-Achse set(gcf,'color','w');% weißer Hintergrund axis([0 fn 0 0.03]);% Skalierung der x- und y-Achse title(['Frequenzbereich (Auflösung: ',num2str(df),' Hz)']);% Titel xlabel('Frequenz in Hz');% Beschriftung x-Achse ylabel('Amplitude in m/s^2');% Beschriftung y-Achse legend('standard Lüfter', 'defekter Lüfter','Location','best','LineWidth',1.0);% Legende grid on;% Gitternetzlinie hold off;% nicht mehr in gleiches Bild plotten % Amplitude in dB % --------------- figure('Name','Bildbereich - Körperschall - standard Lüfter vs. defekter Lüfter - PWM 100 %', ... 'NumberTitle','Off');% Beschriftung Hauptüberschrift plot(x_fs-fn, mag2db((Xfin/bw_a) + eps),'k.-','LineWidth',3.0);% mag2db = 20*log10(xa), eps = kleine Konstante zur Vermeidung von log10(0) hold on;% in gleiches Bild plotten plot(x_fs-fn, mag2db((XMfin/bw_a) + eps),'r.-','LineWidth',1.0);% mag2db = 20*log10(xaM), eps = kleine Konstante zur Vermeidung von log10(0) plot(x_fs-fn, mag2db((Yfin/bw_a) + eps),'k.-','LineWidth',3.0);% mag2db = 20*log10(ya), eps = kleine Konstante zur Vermeidung von log10(0) plot(x_fs-fn, mag2db((YMfin/bw_a) + eps),'r.-','LineWidth',1.0);% mag2db = 20*log10(yaM), eps = kleine Konstante zur Vermeidung von log10(0) plot(x_fs-fn, mag2db((Zfin/bw_a) + eps),'k.-','LineWidth',3.0);% mag2db = 20*log10(za), eps = kleine Konstante zur Vermeidung von log10(0) plot(x_fs-fn, mag2db((ZMfin/bw_a) + eps),'r.-','LineWidth',1.0);% mag2db = 20*log10(zaM), eps = kleine Konstante zur Vermeidung von log10(0) set(gca,'FontSize',12,'XScale','log');% Schriftgröße, logarithmische x-Achse, Skalierung x-Achse set(gcf,'color','w');% weißer Hintergrund axis([0 fn -25 95]);% Skalierung der x- und y-Achse title(['Frequenzbereich (Auflösung: ',num2str(df),' Hz)']);% Titel xlabel('Frequenz in Hz');% Beschriftung x-Achse ylabel('Amplitude in dB');% Beschriftung y-Achse legend('standard Lüfter', 'defekter Lüfter','Location','best','LineWidth',1.0);% Legende grid on;% Gitternetzlinie hold off;% nicht mehr in gleiches Bild plotten %% Mikrofon - Darstellung im Bildbereich % -------------------------------------- % Mikrofon - standard Lüfter % -------------------------- MIC = fft(mic_win, nfft);% FFT mit Fensterfunktion MICabs = abs(MIC);% Absolutzahlen MICfin = fftshift(MICabs/nfft);% Verschiebung der Null-Frequenz-Komponente % zur Mitte des Spektrums % Mikrofon - defekter Lüfter % -------------------------- MICM = fft(micM_win, nfft);% FFT mit Fensterfunktion MICMabs = abs(MICM);% Absolutzahlen MICMfin = fftshift(MICMabs/nfft);% Verschiebung der Null-Frequenz-Komponente % zur Mitte des Spektrums % Graphische Darstellung % ---------------------- % Amplitude in Pa % --------------- figure('Name','Bildbereich - Luftschall - standard Lüfter vs. defekter Lüfter - PWM 100 %', ... 'NumberTitle','Off');% Beschriftung Hauptüberschrift stem(x_fs-fn, MICfin, 'k.-','LineWidth',3.0);% hold on; stem(x_fs-fn, MICMfin, 'r.-','LineWidth',1.0);% set(gca,'FontSize',12,'XScale','log');% Schriftgröße, logarithmische x-Achse, Skalierung x-Achse set(gcf,'color','w');% weißer Hintergrund axis([0 fn 0 3e-5]);% Skalierung der x- und y-Achse title(['Frequenzbereich (Auflösung: ',num2str(df),' Hz)']);% Titel xlabel('Frequenz in Hz');% Beschriftung x-Achse ylabel('Amplitude in Pa');% Beschriftung y-Achse legend('standard Lüfter', 'defekter Lüfter','Location','best','LineWidth',1.0);% Legende grid on;% Gitternetzlinie hold off;% nicht mehr in gleiches Bild plotten % Amplitude in dB % --------------- figure('Name','Bildbereich - Luftschall - standard Lüfter vs. defekter Lüfter - PWM 100 %', ... 'NumberTitle','Off');% Beschriftung Hauptüberschrift plot(x_fs-fn, mag2db((MICfin/bw_p) + eps),'k.-','LineWidth',3.0);% mag2db = 20*log10(x), eps = kleine Konstante zur Vermeidung von log10(0) hold on;% in gleiches Bild plotten plot(x_fs-fn, mag2db((MICMfin/bw_p) + eps),'r.-','LineWidth',1.0);% mag2db = 20*log10(x), eps = kleine Konstante zur Vermeidung von log10(0) set(gca,'FontSize',12,'XScale','log');% Schriftgröße, logarithmische x-Achse, Skalierung x-Achse set(gcf,'color','w');% weißer Hintergrund axis([0 fn -80 35]);% Skalierung der x- und y-Achse title(['Frequenzbereich (Auflösung: ',num2str(df),' Hz)']);% Titel xlabel('Frequenz in Hz');% Beschriftung x-Achse ylabel('Amplitude in dB');% Beschriftung y-Achse legend('standard Lüfter', 'defekter Lüfter','Location','best','LineWidth',1.0);% Legende grid on;% Gitternetzlinie hold off;% nicht mehr in gleiches Bild plotten