function [psd,f,std_f,f_plot]=PSDcalc(Fs,timeSerie,lengthSec)

% PSD calculation   
nfft=round(Fs*lengthSec);
noverlap = round(0.9*nfft);
[psd,f] = psdbk(timeSerie,Fs,nfft,noverlap);
 
% STD calculation   
 I = find(f>f(1) & f<=round(Fs/2.));

tau_ref = length(timeSerie)/Fs;
stdTau_ref=Rmstau(f(I),psd(I),tau_ref);

 % fprintf(' ==> TimeSerie : STD computed with Rmstau = %5.0f unit, std = %5.0f unit \n',stdTau_ref,std(timeSerie));

% STD vs BW calculation
 f_plot=[1:round(Fs/2.)];
 for k=1:length(f_plot)
    I = find(f>f(1) & f<=f(k));
    std_f(k)=Rmstau(f(I),psd(I),tau_ref)/stdTau_ref;
 end

  