function [y1,y2,y3] = banc_de_filtres(x,fe);

% [y1,y2,y3] = banc_de_filtres(x,fe);
% filtrage d'un signal dans les bandes cardiovasculaires VLF, LF, HF
% x signal d'entrée
% fe fréquence d'échantillonnage - paramètre optionnel
% si fe n'est pas rentré, l'affichage est en fréquences normalisées
% y1,y2,y3 : sorties des 3 filtres 

if nargin==1,
    fe = 1;
end
L = length(x);

[N,wn]=cheb1ord(0.04,0.05,0.5,20);
[b1,a1] = cheby1(N,0.5,wn);
[N,wn]=cheb1ord([0.05 0.14],[0.03 0.16],0.5,20);
[b2,a2] = cheby1(N,0.5,wn);
[N,wn]=cheb1ord([0.16 0.39],[0.14 0.41],0.5,20);
[b3,a3] = cheby1(N,0.5,wn);

G1 = freqz(b1,a1,2000);
G2 = freqz(b2,a2,2000);
G3 = freqz(b3,a3,2000);

y1 = filtfilt(b1,a1,x);
y2 = filtfilt(b2,a2,x);
y3 = filtfilt(b3,a3,x);

figure(1)
fr = linspace(0,0.5,2000);
plot(fr*fe,abs(G1),fr*fe,abs(G2),fr*fe,abs(G3))
title('reponse en amplitude des filtres')
if nargin==1,
    xlabel('freq. normalisee')
else 
    xlabel('Hz')
end
subplot(111)

figure(2)
subplot(411)
plot((1:L)/fe,x)
title('signal entree')
if nargin==1,
    xlabel('echantillons')
else 
    xlabel('secondes')
end
subplot(412)
plot((1:L)/fe,y1)
title('sortie filtre VLF')
if nargin==1,
    xlabel('echantillons')
else 
    xlabel('secondes')
end
subplot(413)
plot((1:L)/fe,y2)
title('sortie filtre LF')
if nargin==1,
    xlabel('echantillons')
else 
    xlabel('secondes')
end
subplot(414)
plot((1:L)/fe,y3)
title('sortie filtre HF')
if nargin==1,
    xlabel('echantillons')
else 
    xlabel('secondes')
end
subplot(111)

figure(3)
[Px,f] = pwelch(x-mean(x),L/2,[ ],[ ],fe);
Ind = find(f<=0.4);
subplot(411)
plot(f(Ind),Px(Ind))
title('DSP signal')
[Px,f] = pwelch(y1-mean(y1),L/2,[ ],[ ],fe);
subplot(412)
plot(f(Ind),Px(Ind))
title('DSP VLF')
[Px,f] = pwelch(y2-mean(y2),L/2,[ ],[ ],fe);
subplot(413)
plot(f(Ind),Px(Ind))
title('DSP LF')
[Px,f] = pwelch(y3-mean(y3),L/2,[ ],[ ],fe);
subplot(414)
plot(f(Ind),Px(Ind))
title('DSP HF')
subplot(111)
