function [fmoy,rapport,rapport2] = estim_DSP_3(x,fe,flim);

% [fmoy,rapport,rapport2] = estim_DSP_3(x,fe,flim);
% affiche l'estimation de la DSP du signal x
% fe fréquence d'échantillonnage
% flim limite de fréquence pour la représentation
% sorties:
%   fmoy : frequence moyenne
%   rapport : rapport des puissances (0.05-0.25 Hz)/(0.25-0.5 Hz)
%   rapport2 : rapport puissances pic principal/(0.05-0.25 Hz)

L = length(x);
[Px,f] = pwelch(x-mean(x),round(L/4),[ ],[ ],fe);
Ind = find(f<=flim);

fmoy = f'*Px/sum(Px);

bind = find((f>=0.05)&(f<0.25));
cind = find((f>=0.25)&(f<0.5));
rapport = sum(Px(bind))/sum(Px(cind));
%rapport = sum(Px(bind))/sum(Px);

[bid,imax] = max(Px);
rapport2 = sum(Px(max(1,imax-2):imax+2))/sum(Px(bind));


plot(f(Ind),Px(Ind));
xlabel('Hz')
hold on
vy = get(gca,'Ylim');
plot(fmoy*ones(2,1),vy,'r')
title('bleu: DSP, rouge: freq. moyenne')
hold off




