function [Pxx, f, Pxxc] = psd ( x, Fs, nfft, dflag, window )
% [Pxx, f, Pxxc] = psd ( x, Fs, nfft, dflag, window )
%
% Estimate the positive frequecy part of the double-sided Power Spectral 
% Density of a signal vector X using Welch's averaged periodogram method.  
% The signal X is divided into overlapping segments of length NFFT,
% each of which is detrended, then windowed.  The magnitude squared of the 
% length NFFT DFTs of the segments are averaged to form Pxx.  
% Pxx is length NFFT/2+1 for NFFT even, (NFFT+1)/2 for NFFT odd,
% or NFFT if the signal X is complex. Fs is the sampling frequency which
% doesn't effect the spectrum estimate but is used for scaling of plots.
%
% Pxx is normalized to agree with Parseval's theorem
%        trapz(t,x.^2)/T  equals  2*trapz(f',Pxx)    ...  for 0 < f < Fs/2
% 
%  INPUT   
%   x      vector of sampled  input to the system                  (1 x points)
%   Fs     sampling frequency ... 1 / delta_t                      (1 x 1)
%  nfft    the number of points used in the fast Fourier transform (1 x 1)
%  dflag   detrending type: 'none', 'mean' or 'linear'    default: 'linear'
%  window  window function ... one of ...
%	   'dirichlet', 'tapered', 'sine', 'lanczos', 'hamming', 'hann', 'gauss' 
%          default = 'sine'
%
%  OUTPUT
%   Pxx    the real-valued auto-power spectral density of u      (1 x nfft/2+1)
%   f      the frequency values  corresponding to H              (1 x nfft/2+1)

%       $Revision: 1.12 $  $Date: 1994/01/25 17:59:38 $

x = x(:);		% make sure x is a column vector
n = length(x);		% number of data points

if nargin < 3, nfft = 2^(floor(log2(n))-2);	end

if nargin < 4,	dflag = 'linear'; end	% linear detrending
%		dflag = 'mean';		% average value detrending
%		dflag = 'none';		% no detrending

if nargin < 5,	window = 'sine';	end
% use Sine window for PSD estimates, use Hann window for TFE estimates

if rem(nfft,2) nfft = nfft+1; end   % make nfft an even number

nWindow  = nfft;	% number of points in a segment of data
nOverlap = nWindow / 2;	% number of overlapping points
index    = 1:nWindow;
Pw       = floor(nWindow/10);

if strcmp(window,'dirichlet') window = ones(nWindow,1); end				% Dirichlet 
if strcmp(window,'tapered')   window = [ 0.5*(1-cos(pi*[0:Pw]/Pw)) ones(1,nWindow-2*Pw-2) 0.5*(1+cos(pi*[0:Pw]/Pw)) ]'; end % tapered rect.
if strcmp(window,'sine')      window = sin(pi*(index-0.5)'/nWindow); end		% Sine  
if strcmp(window,'lanczos')   window = sinc(2*(index-1)'/(nWindow-1)-1); end		% Lanczos
if strcmp(window,'hamming')   window = 0.53836 - 0.46164*cos(2*pi*index'/nWindow); end	% Hamming 
if strcmp(window,'hann')      window = (1+cos(2*pi*(index'-nWindow/2)/nWindow))/2; end	% Hann
if strcmp(window,'gauss')     window = exp(-0.5*((index'-0.5-nWindow/2)/(0.2*nWindow)).^2); end % Gauss

K = fix((n-nOverlap)/(nWindow-nOverlap)); % number of windows

if n < nWindow	% zero-pad x if it has length less than the window length
	x(nWindow)=0;  n=nWindow;
end

if 0
	disp(sprintf('   x        = (length %g)',length(x)))
	disp(sprintf('   y        = (length %g)',length(y)))
	disp(sprintf('   nfft     = %g',nfft))
	disp(sprintf('   Fs       = %g',Fs))
	disp(sprintf('   window   = (length %g)',length(window)))
	disp(sprintf('   nOverlap = %g',nOverlap))
	if ~isempty(p)
		disp(sprintf('   p        = %g',p))
	else
		disp('   p        = undefined')
	end
	disp(sprintf('   dflag    = ''%s''',dflag))
	disp('   --------')
	disp(sprintf('   K        = %g',K))
end

% Normalizing scale factor ==> asymptotically unbiased
  KNW2 = K*norm(window)^2;

% Alternative Normalizing scale factor ==> peaks are about right
% KNW2 = K*sum(window)^2;

Pxx = zeros(1,nfft); 

for k=1:K

    if strcmp(dflag,'linear')
        xw = window.*detrend(x(index),1);		% detrend first, then window 
    elseif strcmp(dflag,'mean')
        xw = window.*detrend(x(index),0);
    else 
        xw = window.*(x(index));
    end

    Xx2  = abs(fft(xw)).^2;
    Pxx  = Pxx + Xx2;

    index = index + (nWindow - nOverlap);	% advance to the next segment

end

% Return the positive frequency part of the power spectral density 
if ~any(any(imag(x)~=0)),   % if x is not complex
    if rem(nfft,2),                     % nfft odd
        select = [1:(nfft+1)/2];
    else                                % nfft even
        select = [1:nfft/2+1];		% include DC, include Nyquist
    end
    Pxx = Pxx(select);
else
    select = (1:nfft);
end
f = [select - 1] * Fs / nfft;

% normalize spectral density to agree with Parseval's theorem
% ... trapz(t,x.^2)/T  equals  2*trapz(f,Pxx)

  Pxx = Pxx / ( KNW2 * Fs );

if ( nargout == 3 )            % find confidence interval if needed
    Pxxc = zeros(size(Pxx));

    % Confidence interval from Kay, p. 76, eqn 4.16:
    % (first column is lower edge of conf int., 2nd col is upper edge)
      p = 0.50;               % confidence level
      alfa = 1 - p;
      if (K>0)  RF = 9*K/11; end % reduction factor
      Pxxc = Pxx' * chi2inv([alfa/2 , 1-alfa/2],2*K*RF) / (2*K*RF);
    % confid = 2*K*Pxx*(1./chi2inv([1-alfa/2 alfa/2],2*K)); % wrong. 

end

% ------------------------------------------------------------------ PSD
