
function [y,fre]=psdbk(signal,Fs,nfft,noverlap);

% Computation of PSD (BKO version!) of the variable "signal"
% [y,fre]=psdbk(signal,Fs,nfft,noverlap)
% signal: 	temporal signal (unit)
% Fs: 		Sampling Frequency (Hz)
% nfft: 	number of points in the FFT
% noverlap: 	number of overlaping points in the nfft segments
% fre: 		Frequency vector for the PSD (Hz)
% y: 		PSD of signal (unit^2/Hz)
%
% NOTE: This function (based on MATLAB PSD function), delivers
%       PROPERLY scaled Power Spectral Density in [unit^2/Hz]
%       Single Sided (i.e all energy in the positive frequencies)
%
%  created:	B. Koehler,1992
%  modified	BKO,FEB 97: revised for Matlab 4.2 (use of built-in "psd.m")
%
w=hanning(nfft);		% modify window here, if required
% sc=norm(w)^2/sum(w)^2,
if (nargin < 4), noverlap = 0; end
fre=(Fs/nfft:Fs/nfft:Fs/2)';
% fre=(0:Fs/nfft:Fs/2)';
 p=psd(signal,nfft,Fs,hanning(nfft),noverlap);
 p=p(2:length(p));	% First point is removed (average of the signal)
 y=p*2/Fs;		% scale properly the PSD in [Unit^2/Hz], Single Side (i.e all energy in the positive frequency)
if nfft>length(signal), disp('WARNING!!: nfft more than signal''s length'),disp('Results may be WRONG!!'),end
% loglog(fre,sqrt(y));grid
