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