%
% FILTER_CHIRP_WGN.m
%
% (c) P. Flandrin, 30 septembre 2003 - 23 septembre 2009
%
% calcule, trace et fait entendre la filtrˇe de Wiener 
% d'un chirp linˇaire noyˇ dans un bruit blanc 
%
% calls : - AR_coeff.m 

SNR = 0; % SNR

Fs = 1500;
nb = 1;

% ----- chirp

N = 2^11; % nombre total de points du signal
Ns = N/2;
f1 = .1;%.2; % frˇquence de dˇpart
f2 = .05;%.15; % frˇquence d'arrivˇe
chirp = (hanning(Ns).* real(fmlin(Ns,f1,f2)))';
chirp = chirp/sqrt(sum(chirp.^2));
deb = fix(3*N/8);
chirpeff = zeros(1,N);
chirpeff(deb:deb+Ns-1) = chirp;

figure(nb+1)
clf

subplot(311)

plot(chirpeff);
axis([1 N 1.1*min(chirpeff) 1.1*max(chirpeff)])
set(gca,'YTick',[]);
title('signal')
set(gca,'XTick',[]);

pause

sound(chirpeff/max(abs(chirpeff)),Fs)

pause

% ----- bruit

b = randn(size(chirpeff)) ;
b = b/sqrt(sum(b.^2));

% ----- signal

x = sigmerge(chirpeff',b',SNR);

Te = 1;
t = (0:N-1)*Te;

figure(nb+1)

subplot(312)

plot(t,x)
axis([1 t(N) 1.1*min(x) 1.1*max(x)])
title('signal bruite')
set(gca,'YTick',[]);
set(gca,'XTick',[]);

pause

sound(x/max(abs(x)),Fs)

subplot(312)
plot(t,x)
axis([1 t(N) 1.1*min(x) 1.1*max(x)])
title('signal bruite')
set(gca,'YTick',[]);
set(gca,'XTick',[]);

hold on

[xdec,xfin] = ginput(2);
jnddec = find(t >= xdec(1));
inddec = jnddec(1);
jndfin = find(t >= xdec(2));
indfin = jndfin(1);
ref = x(inddec:indfin);
Nref = length(ref);
plot([xdec(1) xdec(1)],[1.1*min(x) 1.1*max(x)],'r')
plot([xdec(2) xdec(2)],[1.1*min(x) 1.1*max(x)],'r')
hold off

noverlap = 32;
Nfft = 128;
Pb = spectrum(ref,2*N,noverlap,hanning(Nfft))/Nref;

xeff = x(indfin:N);
Nref = length(x);

SNRin = 10*log10(sum(chirpeff.^2)/sum((chirpeff-x').^2));

Pxeff = spectrum(xeff,2*N,noverlap,hanning(Nfft))/Nref;

seuil = .01;
Hf = max(1 - Pb(:,1)./Pxeff(:,1),seuil);

H = zeros(1,2*N);
H(1:N+1) = Hf;
H(N+1:2*N) = Hf(N+1:-1:2);

h = real(ifft(H,2*N));

ind = find(abs(h(1:N))>=h(1)/100);
B = h(ind);
LB = length(B);
BB = [fliplr(B(2:LB)) B];

y = filter(BB,1,x);

Ldec = LB-1;

SNRout = 10*log10(sum(chirpeff.^2)/sum((chirpeff-y').^2));

subplot(313)

plot([y(LB:length(y))' zeros(LB-1,1)']');
axis([1 Nref 1.1*min(y) 1.1*max(y)])
set(gca,'YTick',[]);
title('detail filtre')
set(gca,'XTick',[]);

pause

sound(y/max(abs(y)),Fs)