function [p_ret,f,rms_ret] = psdfft (y,nfft,fsamp,wndw,novlap,max_hold)
% psdfft: Estimates Power Spectral Density by periodogram method using an FFT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Make sure that input is a column vector
if(size(y,1)==1)
   y=y';
end;

%Calculate number of available ensembles
npts = length(y);
ensembles = floor((npts-nfft)/(nfft-novlap))+1;

%Initialize ensemble indexing variables
n1 = 1;
n2 = nfft;
dn = nfft-novlap;

%Create frequency vector
df = fsamp/nfft;
f=[0:df:fsamp/2]';

%Compute window function normalization factor
if (wndw == 1)
   wndw_sc = nfft;
elseif (wndw==2)
   wndw_sc = sum(blackman(nfft).^2);
elseif (wndw==3)
   wndw_sc = sum(hamming(nfft).^2);
elseif (wndw==4)
   wndw_sc = sum(hanning(nfft).^2);
end

%Initialize psd summation storage variable
if (max_hold==0)
   arg_sum=zeros([nfft 1]);
elseif (max_hold==1)
   pmax = zeros([((nfft/2)+1) 1]);
end;

%Main program loop
for k=1:ensembles
   arg_y = y(n1:n2);                           %Extract current ensemble points
   arg_y = arg_y - mean(arg_y);                %Remove mean
   if (wndw == 2)
      arg_y = arg_y .* blackman(nfft)';        %Apply a Blackman window
   elseif (wndw == 3)
      arg_y = arg_y .* hamming(nfft)';         %Apply a Hamming window
   elseif (wndw == 4)
      arg_y = arg_y .* hanning(nfft)';         %Apply a Hanning window
   end
   arg_fft = fft(arg_y,nfft);                  %FFT of ensemble
   arg_abs = abs(arg_fft).^2;                  %Modulus squared of FFT
   if (max_hold==0)
      arg_sum = arg_sum + arg_abs;             %Accumulate in summation variable
   elseif (max_hold==1)
      p = 2*arg_abs(1:nfft/2+1);               %Use lower half, scale
      p = p/wndw_sc;                           %Normalize to correct for window
      p = p/fsamp;                             %Normalize spectral density for sampling units
      pmax = vectormax(pmax, p);               %Find max at each point in vector
   end;
   n1 = n1 + dn;                               %Increment ensemble index
   n2 = n2 + dn;
end

if (max_hold==0)
   %Average value of summed spectra
   arg_sum = arg_sum/ensembles;
   
   %Power spectrum is symmetric about Nyquist frequency, use lower half and multiply by 4.
   % Then divide by 2 to convert peak^2 to rms^2. Thus scale factor is 2.
   p = 2*arg_sum(1:nfft/2+1);
   
   %Normalize to correct for window
   p = p/wndw_sc;
   
   %Normalize spectral density for sampling units
   p = p/fsamp;
   
   %Calculate overall rms level from area under PSD curve
   oarms = sqrt(sum(p.*df));
   
   % values to return
   p_ret=p;
   rms_ret=oarms;
elseif (max_hold==1)
   %Calculate overall rms level from area under PSD curve
   maxrms = sqrt(sum(pmax.*df));
   
   % values to return
   p_ret=pmax;
   rms_ret=maxrms;
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
return;


% |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
function [vmax] = vectormax (currentmax, v)
% vectormax: This function compares each element in "v" with the corresponding element
% in currentmax and updates vmax with the bigger of the two values.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

vmax = currentmax;
for k = 1:length(v)
   if v(k) > currentmax(k)
      vmax(k) = v(k);
   end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
return;


% |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
function [w] = blackman (n)
% blackman: Blackman windowing function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

w = 0.42 - 0.5*cos(2*pi*[0:(n-1)]/(n-1)) + 0.08*cos(4*pi*[0:(n-1)]/(n-1));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
return;


% |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
function [w] = hamming (n)
% hamming: Hamming windowing function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

w = 0.54 - 0.46*cos(2*pi*[1:n]/(n-1));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
return;


% |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
function [w] = hanning (n)
% hanning: Hanning windowing function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

w = 0.5*(1 - cos(2*pi*[1:n]/(n+1)));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
return;
