%%  Nachweis für den Gebrauch der A-Wertung sowie der Umrechnung in dB

%   EDIT!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

%   Zur kompletten Umrechnung der Messdaten, sprich von Rohsignal
%   (Zeitbereich) zum bewerteten Zeitbereich ist es notwendig die richtigen
%   Arbeitsschritte einzuhalten!

%   1.  FFT durchführen
%   2.  A-Bewertung durchführen [(A)]
%   3.  IFFT durchführen
%   4.  Umrechnung zu Schalldruck [dB(A)]

%   Alle anderen Wege haben in einem Versuch keine Vergleichbaren
%   (korrekten) Ergebnisse gezeigt (Vergleichsdaten: XXX)!!!

%%  Einlesen der Daten

%rohdaten1=dlmread('01_white_20kmh_1\23.txt', '\t', 23, 0);

%load Vergleiche_fuer_A_Wertung

%%  Aufteilen der Daten

time1=rohdaten1(:,1);
soundpressurelevel1=(rohdaten1(:,2))';
       
%%  FFT inkl. Betrag, Phase und Skalierung

t11=time1;
x11=soundpressurelevel1';

T11=diff(t11(1:2));
fs11=1/T11;
N11=length(x11);

fn=0.5*fs11;
df=fs11/N11;

f11 = 0:df:fn;

H11=fft(x11);

% Betrag

H11_pos0 = abs(H11(1:(N11/2)+1));

% Phase

phase = angle(H11);

% Skalierung

H11_pos01=[H11_pos0(1)/N11; H11_pos0(2:(N11/2))/(N11/2); H11_pos0(end)/N11];

%%  An FFT die dB Umrechnung vornehmen

%   Es ist darauf zu achten, dass der richtige Rechenweg eingehalten wird
%   um die korrekte Lösung zu erhalten!!!!

%   Die dB-Berechnung im FFT bereich wird nur durchgeführt um die
%   korrektheit der Berechnung mit unserem Vergleichssystem zu
%   kontrollieren und zu beweisen.

%   Um den A-bewerteten Schalldruckpegel zu erhalten ist dieser Schritt
%   nicht notwendig und vorallem FALSCH!!!!!!!!
%   Lp=20*log10(p/po) [dB]  ;    po=0.00002 [Pa]
    
    po=0.00002;

H11_pos01dB=20*log10(H11_pos01/po);

%%  An FFT die A-Bewertung durchgeführt

%   Es ist darauf zu achten, dass erst die FFT und dann die A-Bewertung
%   durchgeführt wird!!!

%   Dieser Schritt ist zur korrekten Berechnung notwendig!!!!!!

%   Ra=((12200^2.*f.^4)./((f.^2+20.6^2).*(f.^2+12200^2).*(sqrt(f.^2+...
%   107.7^2)).*(sqrt(f.^2+737.9^2))));

Ra=((12200^2.*f11.^4)./((f11.^2+20.6^2).*(f11.^2+12200^2).*(sqrt(f11.^2+...
    107.7^2)).*(sqrt(f11.^2+737.9^2))));

%   2 dB Anhebung

Ra=Ra*1.2589254117941673;

H11A=H11_pos01 .* Ra';

%%  An A-Bewertung die dB_umrechnung durchführen

%   Umrechnung (A) in dB(A) im Frequenzbereich ( Lp=20*log10(p/po) [dB]

%   Auch diese Berechnung ist zum Vergleich und zur Fehleranalyse
%   durchgeführt worden.

%   Für Endergebnis nicht notwendig und nicht richtig!!!

H11dBA=20*log10(H11A./po);
H11dBA(1)=0;

%%  Rücktransformation (ifft)

% A-Bewertung ohne Skalierung wird für die Rücktransformation benötigt

H11Ao = H11_pos0.*Ra';

%   beidseitiges Spektrum für ifft wieder zusammensetzen 0...fn [Hz]
%   -fn...-df [Hz]  (Signallänge wieder herstellen)

ReA = [H11Ao;H11Ao(end:-1:2)];

% Umrechnung von Polar- in karthes. Koordinaten

Re = ReA .* cos(phase);
Im = ReA .* sin(phase);

% complex Re + j im

HA_ifft = complex(Re,Im);

% Rücktransformation in den Zeitbereich ( wieder noch in Pa )

soundpressurelevel1_A = ifft(HA_ifft);

%%  dB-Berechnung an rücktransformiertem Signal durchführen

%   Pa in dB

%   Notwendiger Schritt!!!

soundpressurelevel1_dBA=20*log10(soundpressurelevel1_A/po);

%%  Rohdaten in dB( Lp=20*log10(p/po) [dB] // po=0.00002 [Pa]

po=0.00002;
soundpressurelevel1_dB=20*log10(soundpressurelevel1./po);

%%  Plotten aller unter Matlab berechneten Grafiken

%   Es ist mit Sicherheit nicht notwendig alle Grafiken zu plotten!!! Jeder
%   muss für seine Auswertung sich die richtigen Daten ausgeben lassen.
%   Dies ist nur eine Auflistung aller Möglichkeiten, die auch von
%   professioneller Software ausgegeben werden kann!

a=figure('Name','Ablauf der Berechnungen','NumberTitle','off');

    subplot(4,4,1);
        plot(t11,soundpressurelevel1,'r');
        title('Originalsignal');
        xlabel('Zeit [s]'); ylabel('Schalldruck [Pa]');
        axis([0 16 -10 10]); grid on; hold on;
        
    subplot(4,4,2);
        plot(t11,soundpressurelevel1_dB,'r');
        title('Umrechnung Matlab: dB');
        xlabel('Zeit [s]'); ylabel('Schalldruckpegel [Pa(dB)]');
        axis([0 16 -75 150]); grid on; hold on;
        
    subplot(4,4,3);
        plot(t11,soundpressurelevel1_dBA,'r');
        title('Umrechnung Matlab: FFT, (A), ifft, dB(A)');
        xlabel('Zeit [s]'); ylabel('Schalldruckpegel [dB(A)]');
        axis([0 16 -75 150]); grid on; hold on;
        
    subplot(4,4,4);
        plot(f11,H11_pos01dB,'r');
        title('Umrechnung Matlab; FFT, dB');
        xlabel('Frequenz [Hz]'); ylabel('[dB]');
        axis([0 8192 -75 150]); grid on; hold on;
    
    subplot(4,4,5);
        semilogx(f11,H11_pos01,'r');
        title('Umrechnung Matlab: FFT');
        xlabel('Frequenz [Hz]'); ylabel('[/]');
        axis([0 8192 0 0.2]); grid on; hold on;
    
    subplot(4,4,6);
        semilogx(f11,H11A,'r');
        title('Umrechnung Matlab: FFT, (A)');
        axis([0 8192 0 10^-3]); xlabel('Frequenz [Hz]'); ylabel('[(A)]');
        grid on; hold on;
        
    subplot(4,4,7);
        semilogx(f11,H11_pos01dB,'r');
        title('Umrechnung Matlab; FFT, dB');
        xlabel('Frequenz [Hz]'); ylabel('[dB]');
        axis([0 8192 -75 150]); grid on; hold on;
    
    subplot(4,4,8);
        semilogx(f11,H11dBA,'r');
        title('Umrechnung Matlab; FFT, (A), dB(A)');
        xlabel('Frequenz [Hz]'); ylabel('[dB(A)]');
        axis([0 8192 -75 150]); grid on; hold on;
        
    
        
%%  Vergleich Matlab mit LMS
        
    subplot(4,4,9);
        plot(t11,soundpressurelevel1,'r');
        title('Vergleich Originalsignal Matlab => XXX');
        xlabel('Zeit [s]'); ylabel('Schalldruck [Pa]');
        axis([0 16 -10 10]); grid on; hold on;
        
        plot(Rohdaten(:,1),Rohdaten(:,2),'k')
    
    subplot(4,4,10);
        plot(t11,soundpressurelevel1_dB,'r');
        title('Vergleich Umrechnung Matlab: dB => XXX');
        xlabel('Zeit [s]'); ylabel('Schalldruckpegel [Pa(dB)]');
        axis([0 16 -75 150]); grid on; hold on;
        
        plot(RohdatendB(:,1),RohdatendB(:,2),'k')
        
    subplot(4,4,11);
        plot(t11,soundpressurelevel1_dBA,'r');
        title('Vergleich Umrechnung Matlab: FFT, (A), ifft, dB(A) => XXX');
        xlabel('Zeit [s]'); ylabel('Schalldruckpegel [dB(A)]');
        axis([0 16 -75 150]); grid on; hold on;
        
        plot(RohdatendBA(:,1),RohdatendBA(:,2),'k')
        
    subplot(4,4,12);
        plot(f11,H11_pos01dB,'r');
        title('Umrechnung Matlab; FFT, dB => XXX');
        xlabel('Frequenz [Hz]'); ylabel('[dB]');
        axis([0 8192 -75 150]); grid on; hold on;
        
        plot(dBFFTlinear(:,1),dBFFTlinear(:,2),'k');
    
    subplot(4,4,13);
        semilogx(f11,H11_pos01,'r');
        title('Vergleich Umrechnung Matlab: FFT => XXX');
        xlabel('Frequenz [Hz]'); ylabel('[/]');
        axis([0 8192 0 0.2]); grid on; hold on;
        
        semilogx(FFT(:,1),FFT(:,2),'k')
    
    subplot(4,4,14);
        semilogx(f11,H11A,'r');
        title('Vergleich Umrechnung Matlab: FFT, (A) => XXX');
        xlabel('Frequenz [Hz]'); ylabel('[(A)]');
        axis([0 8192 0 10^-3]); grid on; hold on;
        
        semilogx(FFTA(:,1),FFTA(:,2),'k')
        
    subplot(4,4,15);
        semilogx(f11,H11_pos01dB,'r');
        title('Vergleich Umrechnung Matlab; FFT, dB => XXX');
        xlabel('Frequenz [Hz]'); ylabel('[dB]');
        axis([0 8192 -75 150]); grid on; hold on;
    
        semilogx(dBFFT(:,1),dBFFT(:,2),'k');
        
    
    subplot(4,4,16);
        semilogx(f11,H11dBA,'r');
        title('Vergleich Umrechnung Matlab; FFT, (A), dB(A) => XXX');
        xlabel('Frequenz [Hz]'); ylabel('[dB(A)]');
        axis([0 8192 -75 150]); grid on; hold on;
        
        semilogx(dBFFTA(:,1),dBFFTA(:,2),'k')