function [ f1, Sf1, f1_Bin, Sf1_Bin ] = Jason_PSD( Data, df, N, CreatePlot, BinLen )
%function [ f, S1 ] = Jason_PSD( Data, df )
%This function calculates the 1-sided PSD directly using fft
%
%Inputs:
% "Data" is the time series data
% "df" is the frequency step in Hz = 1/tmax
% Optional: "N" is the number of steps to be processed by the fft
% Optional" "CreatePlot" is a flag for plotting
% Optional: "BinLen" is the number of frequency steps to bin the fft by
%
%Outputs:
% "f1" are the 1-sided frequencies in Hz
% "Sf1" is the 1-sided PSD in the units of [Data's unit]^2/Hz
% Optional: "f1_Bin" are the binned 1-sided frequencies in Hz
% Optional: "Sf1_Bin" is the binned 1-sided PSD in the units of [Data's unit]^2/Hz
%
%Internals:
% "NData" is the number of steps in Data
% "N2" is half the number of steps to be processed by the fft
% "domega" is the frequency step in rad/s
% "Data_fft" is the DFT of Data
% "omega1" are the 1-sided frequencies in rad/s
% "Somega1" is the 1-sided PSD in the units of [Data's unit]^2/(rad/s)
% "omega2" are the 2-sided frequencies in rad/s
% "Somega2" is the 2-sided PSD in the units of [Data's unit]^2/(rad/s)
% "f2" are the 2-sided frequencies in Hz
% "Sf2" is the 2-sided PSD in the units of [Data's unit]^2/(Hz)
NData = length(Data);
if nargin < 3 || ( nargin >= 3 && isempty(N) )
N = NData;
end
if mod(N,2) == 1
disp( 'Aborting: N must be even in Jason_PSD.' );
return;
end
N2 = N/2;
domega = 2*pi*df;
Data_fft = fft(Data,N); %truncate Data or pad Data with zeros if N /= length(Data)
%Compute 2-sided psds:
omega2 = domega*[ 0:N2 (-N2+1):(-1) ]';
Somega2 = ( Data_fft/N ).*( conj( Data_fft/N ) )/domega;
f2 = omega2 /( 2*pi );
Sf2 = Somega2*( 2*pi );
%Compute 1-sided psds:
omega1 = omega2 (1:N2);
Somega1 = Somega2(1:N2)*2;
f1 = f2 (1:N2);
Sf1 = Sf2(1:N2)*2;
% Bin the PSDs:
if nargin >= 5 && nargout == 4
MaxLen = length(f1);
Bins = 1:BinLen:MaxLen;
NBins = length(Bins);
f1_Bin = zeros(NBins,1);
Sf1_Bin = zeros(NBins,1);
for iBin = 1:NBins
BinStart = Bins(iBin);
BinIx = BinStart:min(MaxLen,(BinStart+BinLen-1));
f1_Bin (iBin) = mean( f1 (BinIx) );
Sf1_Bin(iBin) = mean( Sf1(BinIx) );
end
elseif nargout ~= 2
disp( 'Aborting: Incorrect arguments in Jason_PSD.' );
return;
end
%Plot the PSDs:
if nargin >= 4 && CreatePlot
figure;
semilogy(f1,Sf1,'b.','DisplayName','Raw');
if nargin >= 5
hold on;
semilogy(f1_Bin,Sf1_Bin,'r-','DisplayName','Binned');
hold off;
end
xlabel('Frequency, Hz')
ylabel('PSD, [Data''s unit]^2/Hz')
legend(gca,'show');
end