clear all;
close all;

fp = 150;
t_1 = 0 : 1/fp : (4 - (1/fp));
t_2 = 4 : 1/fp : (8 - (1/fp));
t_3 = 8 : 1/fp : (12 - (1/fp));
t_4 = 12 : 1/fp : 20;

s_1 = sin(9*t_1*2*pi + 0);
s_2 = sin(10*t_2*2*pi + pi/4); % pi = 180°. 
s_3 = 1*randn(1,length(t_3));
s_4 = sawtooth(t_4*2*pi,0.5);

t = [t_1, t_2, t_3, t_4];
s = [s_1, s_2, s_3, s_4];

figure(1);
subplot(4, 1, 1);
plot(t,s,'r');
title('wygenerowany sygnal');
xlabel('czas [s]');

[xw, fpw] = audioread('mbi04jeden.wav');
xv = decimate(xw, 4);
fpv = fpw / 4;
tv = 0 : 1/fpv : (length(xv)-1)/fpv;

subplot(4, 2, 3);
plot(tv,xv);
title('nagranie');
xlabel('czas [s]');

Nfv = 2 ^ nextpow2( (length(xv)) ); 
frv = fft(xv,Nfv);

fv = linspace(0,fpv,Nfv);
wv = abs(frv);

Nfv21 = Nfv / 2 + 1;
% fv = linspace(0,fpv/2,Nfv21);

% wv_db = 20*log10(wv);

subplot(4,2,4);
plot(fv(1:Nfv21),wv(1:Nfv21));
title('widmo amplitudowe nagrania');
xlabel('czestotliwość [Hz]');

fgr = 1200;
fgrN = fgr / (fpv/2);

M = 31;

h = fir1(M-1, fgrN, 'high');
% h = fir1(M-1, fgrN, 'low');
% h = fir1(M-1, [fgrNd, fgrNu], 'bandpass');
% h = fir1(M-1, [fgrNd, fgrNu], 'stop');

Nfh = 2 ^ nextpow2( M ); 
frh = fft(h,Nfh);

fh = linspace(0,fpv,Nfh);
wh = abs(frh);

Nfh21 = Nfh / 2 + 1;
% fh = linspace(0,fpv/2,Nfh21);

subplot(4,2,5);
plot(fh(1:Nfh21),wh(1:Nfh21));
title('modul transmitancji filtra');
xlabel('czestotliwość unormowana [Hz]');

yv = filter(h,1,xv);

fryv = fft(yv,Nfv);

wyv = abs(fryv);

% wyv_db = 20*log10(wyv);

subplot(4,2,6);
plot(fv(1:Nfv21),wyv(1:Nfv21));
title('widmo amplitudowe po filtracji');
xlabel('czestotliwość [Hz]');

tau_max = 400;
% t_tau_max = 0.1;
% tau_max = round(t_tau_max*fp);
k_s = xcorr(s,s,tau_max);

p_k = -tau_max : 1 : tau_max;
t_k = -tau_max/fp : 1/fp : tau_max/fp;

subplot(4,2,7);
plot(t_k,k_s,'r');
title('autokorelacja');
xlabel('przesuniecie [s]');

k_ss1 = xcorr(s,s_1,tau_max);

subplot(4,2,8);
plot(t_k,k_ss1,'r');
title('korelacja wzajemna');
xlabel('przesuniecie [s]');


t_e0 = 0.15;
pr_e0 = round(t_e0 * fpv);
t_e1 = 0.275;
pr_e1 = round(t_e1 * fpv);

ev = xv(pr_e0:pr_e1);
te = tv(pr_e0:pr_e1)