% close all
clear all                 
% import dat

%%%%%%%%%%%%%%%%%%%%%%%% MERENI HLUKU POZADI %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

load('pozadi','Channel_1_Data')
data=Channel_1_Data;
Fs=32768;
%% 
% vyber zkoumane oblasti signalu

% pocatecni cas intervalu
    t1=1
% koncovy cas intervalu
    t2=4

    Ld=length(data)
    tsu=1/Fs;
    tmaxu=(Ld-1)*tsu;
    tuAxis=0:tsu:tmaxu;
    s1=t1/tmaxu
    s2=t2/tmaxu
    n1=round(s1*(Ld-1))
    n2=round(s2*(Ld-1))

    signal=data(n1:1:n2);
%% 

% prirazeni dat
% y = data;
y = signal;
% delka signalu
L = length(y);  
% definice casove osy
dt=1/Fs;
ts=1/Fs;
tmax=(L-1)*ts;
tAxis=0:ts:tmax;

%% 
% FFT
    nfft = 2^nextpow2(L);
% casove okenko
    y_HannWnd = y.*hanning(L);            
    Ydft_HannWnd = fft(y_HannWnd,nfft)/(L/2);

    mYdft = abs(Ydft_HannWnd);
    mYdft = mYdft (1:nfft/2+1);
    mYdft (2:end-1) = 2* mYdft(2:end-1);
    
%% 
% definice frekvencni osy
f = Fs/2*linspace(0,1,nfft/2+1); 
% prepocet na hladiny
p0=2.1E-5

    for i=1:length(mYdft);
        Lp(i)=20*log10(mYdft(i)/p0);
    end
   
% graf casova a frekvencni domena
  figure(1),
  subplot(2,1,1)
  plot(tAxis,y), grid on
  title('Casova domena pro okoli');
  xlabel('Cas,(s)'); 
  ylabel('p,(Pa)');
  subplot(2,1,2)  
  plot(f,mYdft), grid on 
  axis ([0 3000 0 0.15]); 
  title('Frekvencni spektrum pro okoli');
  xlabel('Frekvence,(Hz)'); 
  ylabel('p,(Pa)');
  
% graf hladina Lp
  figure(2)
  plot(f,Lp), grid on
  axis ([0 3000 0 120]);
  title('Frekvencni spektrum Lp pro okoli');
  xlabel('Frekvence,(Hz)'); 
  ylabel('Lp,(dB)');
  
 
% korekce vahovym filtrem A
 f1 = 20.598997; 
 f2 = 107.65265;
 f3 = 737.86223;
 f4 = 12194.217;
  
f=1000;
Ra=(f4^2*f^4)/((f^2+f1^2)*(f^2+f4^2)*sqrt((f^2+f2^2)*(f^2+f3^2)));
Ra1000=(f4^2*f^4)/((f^2+f1^2)*(f^2+f4^2)*sqrt((f^2+f2^2)*(f^2+f3^2)));
Ka1000=20*log(Ra/Ra1000);

f = Fs/2*linspace(0,1,nfft/2+1);
Ra=(f4^2*f.^4)./((f.^2+f1^2).*(f.^2+f4^2).*sqrt((f.^2+f2^2).*(f.^2+f3^2)));
Ka=20*log10(Ra/Ra1000);

% graf hladiny ak. tlaku vazenou filtrem A
figure(3)
plot(f,Lp+Ka), grid on
  axis ([0 3000 0 120]);
  title('Frekvencni spektrum LpA pro okoli');
  xlabel('Frekvence,(Hz)'); 
  ylabel('LpA,(dB)');




%%%%%%%%%%%%%%%%%%%%%%%% MIKROFON V POZICI 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

load('1','Channel_1_Data')
data=Channel_1_Data;
Fs=32768;
%% 
% vyber zkoumane oblasti signalu

% pocatecni cas intervalu
    t1=1
% koncovy cas intervalu
    t2=4

    Ld=length(data)
    tsu=1/Fs;
    tmaxu=(Ld-1)*tsu;
    tuAxis=0:tsu:tmaxu;
    s1=t1/tmaxu
    s2=t2/tmaxu
    n1=round(s1*(Ld-1))
    n2=round(s2*(Ld-1))

    signal=data(n1:1:n2);
%% 

% prirazeni dat
% y = data;
y = signal;
% delka signalu
L = length(y);  
% definice casove osy
dt=1/Fs;
ts=1/Fs;
tmax=(L-1)*ts;
tAxis=0:ts:tmax;

%% 
% FFT
    nfft = 2^nextpow2(L);
% casove okenko
    y_HannWnd = y.*hanning(L);            
    Ydft_HannWnd = fft(y_HannWnd,nfft)/(L/2);

    mYdft = abs(Ydft_HannWnd);
    mYdft = mYdft (1:nfft/2+1);
    mYdft (2:end-1) = 2* mYdft(2:end-1);
    
%% 
% definice frekvencni osy
f = Fs/2*linspace(0,1,nfft/2+1); 
% prepocet na hladiny
p0=2.1E-5

    for i=1:length(mYdft);
        Lp(i)=20*log10(mYdft(i)/p0);
    end
   
% graf casova a frekvencni domena
  figure(4),
  subplot(2,1,1)
  plot(tAxis,y,'-r'), grid on
  title('Casova domena pro pozici 1');
  xlabel('Cas,(s)'); 
  ylabel('p,(Pa)');
  subplot(2,1,2)  
  plot(f,mYdft,'-r'), grid on 
  axis ([0 3000 0 0.15]); 
  title('Frekvencni spektrum pro pozici 1');
  xlabel('Frekvence,(Hz)'); 
  ylabel('p,(Pa)');
  
% graf hladina Lp
  figure(5)
  plot(f,Lp,'-r'), grid on
  axis ([0 3000 0 120]);
  title('Frekvencni spektrum Lp pro pozici 1');
  xlabel('Frekvence,(Hz)'); 
  ylabel('Lp,(dB)');
  
 
% korekce vahovym filtrem A
 f1 = 20.598997; 
 f2 = 107.65265;
 f3 = 737.86223;
 f4 = 12194.217;
  
f=1000;
Ra=(f4^2*f^4)/((f^2+f1^2)*(f^2+f4^2)*sqrt((f^2+f2^2)*(f^2+f3^2)));
Ra1000=(f4^2*f^4)/((f^2+f1^2)*(f^2+f4^2)*sqrt((f^2+f2^2)*(f^2+f3^2)));
Ka1000=20*log(Ra/Ra1000);

f_1 = Fs/2*linspace(0,1,nfft/2+1);
Ra=(f4^2*f.^4)./((f.^2+f1^2).*(f.^2+f4^2).*sqrt((f.^2+f2^2).*(f.^2+f3^2)));
Ka=20*log10(Ra/Ra1000);

% graf hladiny ak. tlaku vazenou filtrem A
figure(6)
plot(f_1,Lp+Ka,'-r'), grid on
  axis ([0 3000 0 120]);
  title('Frekvencni spektrum LpA pro pozici 1');
  xlabel('Frekvence,(Hz)'); 
  ylabel('LpA,(dB)');

%%%%%%%%%%%%%%%%%%%%%%%% MIKROFON V POZICI 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

load('2','Channel_1_Data')
data=Channel_1_Data;
Fs=32768;
%% 
% vyber zkoumane oblasti signalu

% pocatecni cas intervalu
    t1=1
% koncovy cas intervalu
    t2=4

    Ld=length(data)
    tsu=1/Fs;
    tmaxu=(Ld-1)*tsu;
    tuAxis=0:tsu:tmaxu;
    s1=t1/tmaxu
    s2=t2/tmaxu
    n1=round(s1*(Ld-1))
    n2=round(s2*(Ld-1))

    signal=data(n1:1:n2);
%% 

% prirazeni dat
% y = data;
y = signal;
% delka signalu
L = length(y);  
% definice casove osy
dt=1/Fs;
ts=1/Fs;
tmax=(L-1)*ts;
tAxis=0:ts:tmax;

%% 
% FFT
    nfft = 2^nextpow2(L);
% casove okenko
    y_HannWnd = y.*hanning(L);            
    Ydft_HannWnd = fft(y_HannWnd,nfft)/(L/2);

    mYdft = abs(Ydft_HannWnd);
    mYdft = mYdft (1:nfft/2+1);
    mYdft (2:end-1) = 2* mYdft(2:end-1);
    
%% 
% definice frekvencni osy
f = Fs/2*linspace(0,1,nfft/2+1); 
% prepocet na hladiny
p0=2.1E-5

    for i=1:length(mYdft);
        Lp_2(i)=20*log10(mYdft(i)/p0);
    end
   
% graf casova a frekvencni domena
  figure(7)
  subplot(2,1,1)
  plot(tAxis,y), grid on
  title('Casova domena pro pozici 2');
  xlabel('Cas,(s)'); 
  ylabel('p,(Pa)');
  subplot(2,1,2)  
  plot(f,mYdft), grid on 
  axis ([0 3000 0 0.15]); 
  title('Frekvencni spektrum pro pozici 2');
  xlabel('Frekvence,(Hz)'); 
  ylabel('p,(Pa)');
  
% graf hladina Lp
  figure(8)
  plot(f,Lp_2), grid on
  axis ([0 3000 0 120]);
  title('Frekvencni spektrum Lp pro pozici 2');
  xlabel('Frekvence,(Hz)'); 
  ylabel('Lp,(dB)');
  
 
% korekce vahovym filtrem A
 f1 = 20.598997; 
 f2 = 107.65265;
 f3 = 737.86223;
 f4 = 12194.217;
  
f=1000;
Ra=(f4^2*f^4)/((f^2+f1^2)*(f^2+f4^2)*sqrt((f^2+f2^2)*(f^2+f3^2)));
Ra1000=(f4^2*f^4)/((f^2+f1^2)*(f^2+f4^2)*sqrt((f^2+f2^2)*(f^2+f3^2)));
Ka1000=20*log(Ra/Ra1000);

f_2 = Fs/2*linspace(0,1,nfft/2+1);
Ra=(f4^2*f.^4)./((f.^2+f1^2).*(f.^2+f4^2).*sqrt((f.^2+f2^2).*(f.^2+f3^2)));
Ka_2=20*log10(Ra/Ra1000);

% graf hladiny ak. tlaku vazenou filtrem A
figure(9)
plot(f_2,Lp_2+Ka_2), grid on
  axis ([0 3000 0 120]);
  title('Frekvencni spektrum LpA pro pozici 2');
  xlabel('Frekvence,(Hz)'); 
  ylabel('LpA,(dB)');
  
%%%%%%%%%%%%%%%%%%%%%%%% POROVNANI JEDNOTLIVYCH MERENI %%%%%%%%%%%%%%%%%%%%

figure(10)
plot(f_1,Lp+Ka,'-r',f_2,Lp_2+Ka_2,'-b'), grid on
legend('Pozice mikrofonu 1','Pozice mikrofonu 2',5)
axis ([0 3000 0 120]);
  title('Frekvencni spektrum LpA. Srovnani pozic 1 a 2.');
  xlabel('Frekvence,(Hz)'); 
  ylabel('LpA,(dB)');
