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