%
% FA_CHIRP.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

SNR = 0 ; % SNR 50, 5, 0, -5, -10

figure(1)
clf 

hd = 0; hf = 0; % gabarit

Ns = 200; % 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)));

sig = sigmerge(reff,randn(size(reff)),SNR-20*log10(3));%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(3,1,1)
plot(real(reff),'k')
axis([1+Ns dNs+Ns -ampp ampp])
set(gca,'YTick',[]);set(gca,'XTick',[])
hold on
plot([1+Ns dNs+Ns],[0 0],'b')
plot([Ns+Ns+1 Ns+Ns+1],[-ampp ampp],':')
title(['noise-free signal']);%,'FontSize',15)
% xlabel('time','FontSize',15)
hold off
pause

subplot(3,1,2)
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],'b')
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(3,1,2)
	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(3,1,3)
	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(.01)

end