%
% FA_BAT.m
%
% (c) P. Flandrin, 23 janvier 2002 - 24 septembre 2003
%
% calcule et trace en ligne la sortie (enveloppe au carrŽ)
% du filtre adaptŽ ˆ un chirp d'Žcholocation

figure(1)
clf 

Fs = 4196;

load veille.mat
Te = 1/230.4; % ms

xx = veille(1:length(veille));
[B,A] = butter(4,.8);
x = filter(B,A,xx);
N = length(x);
t = (0:N-1)*Te;

subplot(211)

plot(t,x)
axis([1 t(N) 1.1*min(x) 1.1*max(x)])
xlabel('temps (ms)')
title('signal d echolocation + reverberation')
set(gca,'YTick',[]);
hold on

pause

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

[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

pause

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

pause

y = filter(fliplr(ref),1,x);

subplot(212)

ytr = sqrt(abs(y));

plot(t,ytr);
axis([1 t(N) 0 1.1*max(ytr)])
set(gca,'YTick',[]);
title('filtre adapte')
xlabel('temps (ms)')

zoom on

% hd = 0; hf = 0; % gabarit
% 
% Ns = 150; % nb de points du signal utile
% f1 = .5; % frŽquence de dŽpart du chirp
% f2 = .0; % frŽquence d'arrivŽe du chirp
% ref = fmlin(Ns,f1,f2).*amgauss(Ns,Ns/2,Ns/4);
% reff = [zeros(size(ref))' ref' zeros(size(ref))']'; % signal zero-padded
% %nor = abs(sum(reff.*conj(reff)));
% 
% SNR = 30; % SNR
% sig = sigmerge(reff,randn(size(reff)),SNR);%reff + hilbert(randn(size(reff)))/2.5; % signal complet bruitŽ
% amp = 1.2*max(max(abs(sig)),abs(min(abs(sig))));
% ampp = amp*2;
% 
% dNs = 3*Ns/2-1;%2*Ns-1; % horizon temporel de visalisation
% 
% subplot(2,1,1)
% plot(real(sig))
% axis([1+Ns dNs+Ns -ampp ampp])
% set(gca,'YTick',[]);set(gca,'XTick',[])
% hold on
% plot([1+Ns dNs+Ns],[0 0],'r')
% plot([Ns+Ns+1 Ns+Ns+1],[-ampp ampp],':')
% title(['signal + noise (',int2str(SNR),' dB)'],'FontSize',15)
% xlabel('time','FontSize',15)
% hold off
% pause
% 
% cor = [];
% w = [];
% 
% f = 1.5;
% 
% for k = Ns+1:Ns+1+dNs
% 
% 	w = zeros(3*Ns,1);
% 	w(k-Ns:k-1) = ref;
% 	cor = [cor abs(sum(sig.*conj(w)))];
% 
% end
% 
% nor = max(abs(cor));
% 
% cor = [];
% w = [];
% 		
% for k = Ns+1:Ns+1+dNs
% 
% 	subplot(2,1,1)
% 	plot(real(sig),'--')
% 	title(['signal + noise (',int2str(SNR),' dB)'],'FontSize',15)
% 	amp = 1.2*max(max(real(sig)),abs(min(real(sig))));
% 	axis([1+Ns dNs+Ns+1 -ampp ampp])
% 	set(gca,'YTick',[]);set(gca,'XTick',[])
% 	hold on
% 	plot([Ns+Ns+1 Ns+Ns+1],[-ampp ampp],':')
% 	plot([k k],[-ampp ampp],'k')
% 	w = zeros(3*Ns,1);
% 	w(k-Ns:k-1) = ref;
% 	plot(real(w),'k')
% 	hold off
% 		
% 	subplot(2,1,2)
% 	cor = [cor abs(sum(sig.*conj(w)))];
% 	plot(cor.^2,'r')
% 	hold on
% 	plot(-cor.^2,'r')
% 	plot([k-Ns k-Ns],[-f*nor.^2 f*nor.^2],'r')
% 	plot([Ns+1 Ns+1],[-f*nor.^2 f*nor.^2],':')
% 	plot([1 dNs],[0 0])
% 	hold off
% 	axis([1 dNs -1.2*nor.^2 1.2*nor.^2])
% 	set(gca,'YTick',[]);set(gca,'XTick',[])
% 	title('quadrature matched filter output','FontSize',15)
% 	xlabel('time','FontSize',15)
%         
% 	pause(.1)
% 
% end
%