% ANNUL_ECHO.M
%
% P. Flandrin, Sept. 28, 2003 --- Feb. 13, 2006
%
% calcule et trace un signal avec 残ho (+ bruit),
% estime le d残alage et l'att始uation par auto-corr四ation
% et reconstruit par filtrage inverse

% ----- signal

% load handel.mat; x = y; Fs = 8192; dec = 6000; % 6000
% load gong.mat; x = y; Fs = 8192; dec = 1200; % 1200
% load laughter.mat; x = y; Fs = 8192; dec = 1750 ;
% load splat.mat; x = y; Fs = 8192; dec = 750; % 750
% load train.mat; x = y; Fs = 8192; % contre-exemple bande 師roite
% load chirp.mat; x = y; Fs = 8192; dec = 750; % 750
% anna = wavread('annak_soleil'); x = anna(300001:400000); Fs = 22300; dec = 8000; % 5000
 zoe = auread('zoe'); x = zoe; Fs = 44600; dec = 8000 ;

N = length(x);

% N = 2000; x = real(hanning(N).*fmlin(N,0,.5)); Fs = 2048; dec = 15; % 15

figure(1)
clf

subplot(221)
plot(x)
axis([1 N -1.1*max(abs(x)) 1.1*max(abs(x))])
title('signal original')
set(gca,'XTick',[])
set(gca,'YTick',[])

pause

% sound(x,Fs);
p = audioplayer(x, Fs); 
play(p);

pause

% ----- 残ho

a = 0.8 ; % att始uation de l'残ho
SNR = 200;% SNR en dB

B = zeros(1,dec);
B(1) = 1;
B(dec) = a;

yy = filter(B,1,x);
y = sigmerge(yy,randn(size(yy)),SNR);

subplot(222)
plot(y)
axis([1 N -1.1*max(abs(y)) 1.1*max(abs(y))])
hold on
plot([dec dec],[-1.1*max(abs(y)) 1.1*max(abs(y))],'r')
title('signal avec echo')
set(gca,'XTick',[])
set(gca,'YTick',[])
hold off

pause

% sound(y,Fs);
p = audioplayer(y, Fs); 
play(p);

pause

% ----- estimation du d残alage et de l'att始uation

c = xcorr(y);

subplot(224)
plot(N:N+2*dec,c(N:N+2*dec))
axis([N N+2*dec 1.1*min(c(N:N+2*dec)) 1.1*max(c(N:N+2*dec))])
set(gca,'XTick',[])
set(gca,'YTick',[])
title('auto-correlation')
hold on
plot([N+dec N+dec],[.8*1.1*max(c(N:N+2*dec)) 1.1*max(c(N:N+2*dec))],'r')

[xdec1,ydec1] = ginput(2);
xdec = fix(xdec1);
ind = find(abs(c(xdec(1):xdec(2))) == max(abs(c(xdec(1):xdec(2)))));
decest = xdec(1)-N+ind(1);
r = c(N)/c(xdec(1)+ind(1)-1);
% r = cmaxsort(1)/cmaxsort(2);
aest = real(.5*(r-sign(r)*sqrt(r^2-4)));
[decest aest]

text(fix(N+dec/4),.8*1.1*max(c(N:N+2*dec)),['dec = ',num2str(dec)])
text(fix(N+dec/4),.6*1.1*max(c(N:N+2*dec)),['est = ',num2str(decest)])
text(fix(N+dec+dec/4),.8*1.1*max(c(N:N+2*dec)),['att = ',num2str(a,2)])
text(fix(N+dec+dec/4),.6*1.1*max(c(N:N+2*dec)),['est = ',num2str(aest,2)])

hold off

% decest = dec-1;
% aest = .7447;

% ----- filtrage inverse

Best = zeros(1,decest);
Best(1) = 1;
Best(decest) = aest;

xest = filter(1,Best,y);

subplot(223)
plot(xest)
axis([1 N -1.1*max(abs(xest)) 1.1*max(abs(xest))])
title('signal reconstruit')
hold on
plot([dec dec],[-1.1*max(abs(xest)) 1.1*max(abs(xest))],'r')
plot(x-xest,'r')
set(gca,'XTick',[])
set(gca,'YTick',[])
hold off

pause

% sound(xest,Fs)
p = audioplayer(xest, Fs); 
play(p);