clear variables;
close all;

[s,Fs] = audioread('ProtestMonoBruit.wav');

figure(1);
N = length(s);
Ts = 1/Fs;
temps = 0:Ts:(N-1)*Ts;
plot(temps,s);

% sound(s,Fs);

figure(2);
temps2 = temps(60*Fs:70*Fs);
s2 = s(60*Fs:70*Fs);
[R,lags] = xcorr(s2,200,'biased');
plot(lags,R);

figure(3);
M = 20;
origine = ceil(length(R)/2);
B = zeros(1,M+1);
B(1) = 1;
V_gamma = R(origine:origine+M);
M_gamma = toeplitz(V_gamma);
Minv_gamma = pinv(M_gamma);
Phi = Minv_gamma*B';
sigma = sqrt(1/Phi(1));
H = -sigma^2*Phi;
H(1) = 0;
bar(H);

figure(4);
s2_estime = [];
for n = M+1:100000;
    s2_lags = s2(n:-1:n-M);
    s2_estime = [s2_estime,H'*s2_lags];
end
subplot(211);
temps2_decale = temps2(M+1:100000);
s2_decale = s2(M+1:100000)';
plot(temps2_decale,s2_estime,'r',temps2_decale,s2_decale);


subplot(212);
erreur = abs(s2_decale-s2_estime);
plot(temps2_decale,erreur);

figure(5);
seuil = 0.03;
s2_bis = s2;
for l = find(abs(erreur) > seuil)
    s2_bis(l) = median(s2(l-10:l+10));
end
plot(temps2,s2_bis,'r',temps2,s2);

figure(6);
subplot(211);
plot(temps2,s2);
subplot(212);
erreur2 = abs(s2_bis-s2);
plot(temps2,erreur2);
axis([60 63 0 0.1]);