function [f,psd,ENBW,DC] = mypsd(x,fs,N,window,overlap)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% [f,psd,ENBW] = mypsd(x,fs,N,window,overlap)
% 
% Computes Power Spectral Density using Welch's method of 
% overlapped segmented averaging of modified periodograms.
% Based on techniques described in "Spectrum and spectral density
% estimation by the Discrete Fourier Transform", G. Heinzel, et. al.
% 
% INPUTS:
%   x = time series input [amplitude]
%   fs = sampling frequency [Hz]
%   N = desired length of FFT 
%      NOTE: bin spacing = lowest frequency = fs/N
%   window = type of window ('flat', 'hanning', etc.)
%      Suported Windows:
%         - none (rectangular, boxcar)
%         - Bartlett (triangle)
%         - Hanning
%         - Hamming
%         - Blackmann-Harris (92dB)
%         - Kaiser (kaiser3, kaiser4, kaiser5)
%   overlap = fractional overlap between sections 
%      NOTE: use overlap = 'auto' to use reccommended 
%            overlap for chosen window.
%
% OUTPUTS:
%   f = vector of Fourier frequencies
%   psd = power spectral density [amplitude^2/Hz]
%   ENBW = equivalent noise bandwidth [Hz]
%      NOTE: Other useful spectra can be determined as follows
%         - Power Spectrum, ps= psd * ENBW [amplitude^2]
%         - Linear Spectral Density, lsd = sqrt(psd) [amplitude/sqrt(Hz)]
%         - Linear Spectrum, ls = sqrt(ps) [amplitude]
%
% Ira Thorpe
% Updated 4-26-05
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%disp(' ');

j = 0:N-1;                % points for time series
f_res = fs/N;             % frequency resolution
t = (1/fs).*j;            % time values (used in linear regression)
f = (0:N/2).*f_res;       % vector of Fourier frequencies

% Definition of windows (from G. Heinzel, et. al.)
switch lower(window)
    case {'rectangular','flat','boxcar','none'}
        w = ones(1,N);
        ROV = 0;
        %NENBW = 1;
        disp('No Windowing');
    case {'bartlett','triangle'}
        z = 2.*j/N;
        for i = 1:length(z)
            if z(i) <= 1
                w(i) = z(i);
            else
                w(i) = 2-z(i);
            end
        end
        ROV = 0.5;
        %NENBW = 4/3;
        disp('Window is Bartlett');
    case 'hanning'
        z = 2*pi/N.*j;
        w = (cos((z-pi)/2).^2);
        ROV = 0.5;
        %NENBW = 1.5;
        %disp('Window is Hanning');
    case 'hamming'
        z = 2*pi/N.*j;
        w = 0.54-0.46.*cos(z);
        ROV = 0.5; 
        %NENBW = 1.3628;
        disp('Window is Hamming');
    case {'blackman','blackman-harris'}
        % minimum 4-term B-H window
        % or 92dB B-H window
        z = 2*pi/N.*j;
        w = 0.35875-0.48829.*cos(z)+0.14128.*cos(2.*z)-0.01168.*cos(3.*z);
        ROV = 0.661;
        %NENBW = 2.0044;
        disp('Window is Blackman-Harris');
    case 'hft116d'
        z = 2*pi/N.*j;
        w = 1-1.9575375*cos(z)+1.4780705*cos(2*z)-0.6367431*cos(3*z)+0.1228389*cos(4*z)-0.0066288*cos(5*z);
        %NENBW = 4.2186;
        ROV = .782;
        disp('Window is HFT116D');
    case 'kaiser3'
        z = 2/N.*j-1;
        alpha = 3;
        w = besseli(0,pi*alpha.*sqrt(1-z.^2))/besseli(0,pi*alpha);
        %NENBW = 1.7952;
        ROV = 0.619;
        disp('Window is Kaiser, alpha = 3');
    case 'kaiser4'
        z = 2/N.*j-1;
        alpha = 4;
        w = besseli(0,pi*alpha.*sqrt(1-z.^2))/besseli(0,pi*alpha);
        %NENBW = 2.0533;
        ROV = 0.67;
        disp('Window is Kaiser, alpha = 4');
    case 'kaiser5'
        z = 2/N.*j-1;
        alpha = 5;
        w = besseli(0,pi*alpha.*sqrt(1-z.^2))/besseli(0,pi*alpha);
        %NENBW = 2.2830;
        ROV = 0.705;
        disp('Window is Kaiser, alpha = 5');
        
        
    otherwise
        error('Unknown window fucntion');
        
end

% compute window sums and ENBW
S1 = sum(w);
S2 = sum(w.^2);
NENBW = N*S2/(S1^2);
ENBW = NENBW*f_res; 
 
% Determine overlap
if lower(overlap) == 'auto'
    overlap = ROV; 
elseif overlap > 1 || overlap < 0
    error('Overlap must be between 0 and 1');
end
%disp(['Overlap is ', num2str(overlap*100), '%']);

% Determine PSD 
%disp('Computing PSD...');
tic 

Navg = floor(length(x)/N*(1+overlap));      % # of segments
% disp(['Number of averages = ', num2str(Navg)]);
 
for i = 1:Navg;   
    j_low = floor(N*(i-1)*(1-overlap)+1);                % lower index for segment i
    j_high = j_low+N-1;                                  % upper index for segment i
    seg = x(j_low:j_high);                               % segment i
    fit = polyfit(t,seg,1);                              % linear fit to segment i
    DC = polyval(fit,t);
    res = seg - DC;                          % residual from linear fit for segment i
    win_res = res.*w;                                    % windowed residual
    y = fft(win_res);                                    % DFT of windowed residual
    p(i,:) = 2*y(1:N/2+1).*conj(y(1:N/2+1))/(fs*S2);     % PSD for segment i

end
psd = mean(p,1);    % Average over all segments to get PSD for entire data set
%disp(['Finished. Elapsed time = ', num2str(toc),' seconds.']);
 
% END PROGRAM