% detect_tf.m
%
% P. Flandrin
clf 

distr = 1; % choix de la distr. tf: PWV = 1; SP = 2; RSP = 3
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

sig = 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(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','FontSize',16)
xlabel('time','FontSize',16)
hold off
pause

if distr == 1

	wv = tfrpwv(sig,Ns+1:1+dNs+Ns,128,window(127,'Kaiser'),1);
	tf = wv+abs(wv);
	subplot(312)
	title('Wigner-Ville','FontSize',16)

elseif distr == 2

	sp = tfrsp(sig,Ns+1:1+dNs+Ns,256,window(31,'Kaiser'),1);
	sp = sp/max(max(sp));
	tf = sp;
	subplot(312)
	title('spectrogram','FontSize',16)

else

	[sp,rsp] = tfrrsp(sig,Ns+1:1+dNs+Ns,256,window(31,'Kaiser'),1);
	rsp = rsp/max(max(rsp));
	tf = rsp;
	subplot(312)
	title('reassigned spectrogram','FontSize',16)

end

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,1)
	plot(real(sig),'--')
	title('signal + noise','FontSize',16)
	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,2)
	imagesc(flipud(tf(1:128,:)))
	hold on
	plot([k-Ns+1 k-Ns+1],[1 128],'w')
	plot([k-Ns+1-Ns/2-Ns/4+Ns/8 k-Ns+1-Ns/2+Ns/4-Ns/8],[64-32+16-hf 64+32-16-hd],'w','LineWidth',4)
	plot([k-Ns+1-Ns/2-Ns/4+Ns/8 k-Ns+1-Ns/2+Ns/4-Ns/8],[64-32+16+hf 64+32-16+hd],'w','LineWidth',4)
	plot([k-Ns+1-Ns/2-Ns/4+Ns/8 k-Ns+1-Ns/2-Ns/4+Ns/8],[64-32+16-hf 64-32+16+hf],'w','LineWidth',4)
	plot([k-Ns+1-Ns/2+Ns/4-Ns/8 k-Ns+1-Ns/2+Ns/4-Ns/8],[64+32-16-hd 64+32-16+hd],'w','LineWidth',4)
	hold off
	set(gca,'YTick',[]);set(gca,'XTick',[])
	ylabel('frequency','FontSize',16)
		
	if distr == 1

		title('Wigner-Ville','FontSize',16)

	elseif distr == 2

		title('spectrogram','FontSize',16)

	else

		title('reassigned spectrogram','FontSize',16)

	end

	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',16)
	xlabel('time','FontSize',16)
        
	pause(.1)

end

