% faTFdopp.m
%
% P. Flandrin
Ns = 200; % nb de points du signal utile
f1 = .05; % frŽquence de dŽpart du chirp
f2 = .45; % frŽquence d'arrivŽe du chirp
a = .7; % facteur Doppler

refl = fmlin(Ns,f2,f1).*amgauss(Ns,Ns/2,Ns/4);
refl = refl/max(abs(refl).^1);
refla = fmlindopp(Ns,f2,f1,Ns/2,a).*amgaussdopp(Ns,Ns/2,Ns/4,a);
refla = refla/max(abs(refla).^1);
reffl = [zeros(size(refl))' refla' zeros(size(refl))']'; % signal zero-padded

alpha = 300;
refh = hilbert(altes(Ns,f1,f2,alpha));
refh = refh/max(abs(refh).^1);
refha = hilbert(altesdopp(Ns,f1,f2,alpha,a));
refha = refha/max(abs(refha).^1);
reffh = [zeros(size(refh))' refha' zeros(size(refh))']'; % signal zero-padded

noise = hilbert(randn(size(reffl)))/2;

figure(1)
clf 

sig = reffl + noise; % signal complet bruitŽ
amp = 1.2*max(max(abs(sig)),abs(min(abs(sig))));
ampp = amp;

dNs = 3*Ns/2-1; % horizon temporel de visalisation

subplot(4,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('Dopplerized linear chirp','FontSize',16)
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) = refl;
	cor = [cor abs(sum(sig.*conj(w)))];

end

nor = max(abs(cor));

cor = [];
w = [];
		
for k = Ns+1:Ns+1+dNs

    subplot(4,1,1)
    plot(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) = refl;
    plot(real(w),'k')
    title('Dopplerized linear chirp','FontSize',20)
    hold off

    subplot(4,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',16)
        
    pause(.1)

	end

pause

sig = reffh + noise; % signal complet bruitŽ
amp = 1.2*max(max(abs(sig)),abs(min(abs(sig))));
ampp = amp;

subplot(4,1,3)
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('Dopplerized hyperbolic chirp','FontSize',16)
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) = refh;
	cor = [cor abs(sum(sig.*conj(w)))];

end

nor = max(abs(cor));

cor = [];
w = [];
		
for k = Ns+1:Ns+1+dNs

    subplot(4,1,3)
    plot(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) = refh;
    plot(real(w),'k')
    title('Dopplerized hyperbolic chirp','FontSize',16)
    hold off

    subplot(4,1,4)
	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',16)
        
    pause(.1)

end

