function [f,PSD,OASPL] = psdcalcfast(x,fs,ns,N)

% This program calulates the FFT, PSD, and OASPL of a signal. This version
% is the exact same as psdcalc2, but instead of using a for loop to make
% your matrix of different blocks, it uses indexing. This is a huge
% improvement, time-wise, if N>>ns. However, this version takes out the
% possibility of using pad.
%
%   call [f,PSD,OASPL] = psdcalc(x,fs,ns,N);
% 
%   Outputs: 
%   f = frequency array for plotting Xss.
%   PSD = Power spectral density
%   OASPL = Overall sound pressure level
%
%   Inputs:
%   x = time series data.
%   fs = sampling frequency
%   ns = number of samples per block.  Default is 2^15 if not specified.
%   N = total number of samples.  If N is not an integer multiple of ns, 
%       the samples less than ns in the last block are discarded.  Default   
%       is nearest lower power of 2 if not specified.
%
%   Authors: Kent Gee and Alan Wall


warning off

%DEFAULT PAD
pad=1;
%{
if nargin<5
    pad=1;
end
%}

%DEFAULT WAVEFORM SIZE
if nargin<4
    N = 2^floor(log2(length(x))); 
end
x = x(1:N);

%DEFAULT BLOCK SIZE
if nargin<3
    ns = 2^15; 
end

% Data pad

padlength=(pad-1)*ns;
padx=zeros(padlength,1);
nspad=ns+padlength;


%FREQUENCY ARRAY
f = fs*(0:nspad/2-1)/nspad;
df = f(3) - f(2);   %Width of frequency bins.



%Enforce zero-mean
x = x-mean(x);


%HANNING WINDOW
ww = hann(ns);
W = mean(ww.*conj(ww)); %Used to scale the psd



%SPLITS DATA INTO BLOCKS
% Divides total data set into blocks of length ns with 50% overlap, and
% windowed.
blocklen = floor(2*N/ns-1);

blockmat = repmat(1:ns,blocklen,1) + repmat(ns/2*(0:blocklen-1)',1,ns);
blocks = repmat(ww',blocklen,1).*x(blockmat);                                               

%COMPLEX PRESSURE AMPLITUDE
X = sqrt(2*df/ns/fs/W)*fft(blocks,nspad,2);    %Correct values for a single-
                                            %sided fft.   
Xss = X(:,1:nspad/2);  %Takes first ns/2 points to make it single-sided. 
                    %Units are Pa.

%POWER SPECTRAL DENSITY
    PSD =(1/df)*mean(Xss.*conj(Xss),1); %Units are Pa^2/Hz

%OVERALL SOUND PRESSURE LEVEL
if nargout > 2
    OASPL = 20*log10(sqrt(sum(PSD*df))/2e-5); 
end

end
