% FA_VIRGO.M
%
% P. Flandrin, Oct. 5, 2003
%
% calcule et trace la sortie du filtre adaptŽ ˆ
% une onde type Virgo
%
% calls : - bhchirp.m

% ----- chirp Virgo

m=30; 
r=100; 
tin= 3; 
fs=5e3;
[t,x,f]=bhchirp(m,r,tin,fs);
x = x/sqrt(sum(x.^2));
N = length(x);

figure(1)
clf

subplot(311)
plot(x)
hold on
plot([1 N],[0 0])
axis([1 N 1.1*min(x) 1.1*max(x)])
set(gca,'XTick',[])
set(gca,'YTick',[])
title('chirp de binaire coalescente')

pause

% ----- rŽfŽrence pour le filtre adaptŽ

ind = find(x == max(x));
Lh = 2500; % longueur du filtre (# Žchantillons avant coalescence)
h = x(ind-Lh:ind);

plot(ind-Lh:ind,x(ind-Lh:ind),'r')
plot([ind-Lh ind-Lh],[1.1*min(x) 1.1*max(x)],'r')
plot([ind ind],[1.1*min(x) 1.1*max(x)],'r')
title('chirp de binaire coalescente + reference pour le filtre adapte')

hold off

pause

% ----- observation

SNR = -10; % SNR
% xb = sigmerge(x,randn(size(x)),SNR);
b = randn(size(x));
b = b/sqrt(sum(b.^2));
alpha = 10^(SNR/20);
xb = b + x*alpha;

subplot(312)
plot(xb)
axis([1 N 1.1*min(xb) 1.1*max(xb)])
set(gca,'XTick',[])
set(gca,'YTick',[])
title(['observation bruitee, SNR = ',int2str(SNR),' dB'])
Ax = axis;

pause

% ----- filtrage adaptŽ

y = filter(flipud(h),1,xb);

btest = randn(size(xb));
btest = btest/sqrt(sum(btest.^2))*sqrt(sum(xb.^2));
yb = filter(flipud(h),1,btest);
ybs = sort(abs(hilbert(yb)));
seuil95 = ybs(fix(length(yb)*.95));
seuil99 = ybs(fix(length(yb)*.999));

subplot(313)
plot(abs(hilbert(y)),'r')
jnd = find(abs(hilbert(y)) == max(abs(hilbert(y))));
hold on
plot([1 N],[seuil95 seuil95],'k')
plot([1 N],[seuil99 seuil99],'k--')
plot([jnd jnd],[0 1.1*max(y)],'r')
plot([ind ind],[0 1.1*max(y)],'r--')
axis([1 N 0 1.1*max(y)])
set(gca,'XTick',[])
set(gca,'YTick',[])
title('enveloppe de la sortie du filtre adapte')
hold off