function y=rmspsd(fre,S,fmin,fmax)
% y=rmspsd(fre,S,fmin,fmax) computes the rms (square root of the power)
% associated with the Power Spectral Density S, between the frequency
% limits fmin and fmax.
%
% f: frequency vector 		(Hz)
% S: PSD			(unit^2/Hz)
% fmin: lower frequency limit 	(Hz)
% fmax: upper frequency limit 	(Hz)
% y : rms			(unit) 
%
ii=find(fmin<=fre & fre<=fmax);
imin=ii(1);imax=ii(length(ii));
df=fre(imin+1:imax)-fre(imin:imax-1);
ds=(S(imin:imax-1)+S(imin+1:imax))/2;
var=sum(ds.*df);
y=sqrt(var);

