function [f, Spec] = jfft(x, nfft, Fs, window, noverlap, LOG)
%PSD Power Spectral Density estimate.
% Pxx = PSD(X,NFFT,Fs,WINDOW) estimates the Power Spectral Density of
% signal vector X using Welch's averaged periodogram method. X is
% divided into overlapping sections, each of which is detrended, then
% windowed by the WINDOW parameter, then zero-padded to length NFFT.
% The magnitude squared of the length NFFT DFTs of the sections are
% averaged to form Pxx. Pxx is length NFFT/2+1 for NFFT even, (NFFT+1)/2
% for NFFT odd, or NFFT if the signal X is complex. If you specify a
% scalar for WINDOW, a Hanning window of that length is used. Fs is the
% sampling frequency which doesn't affect the spectrum estimate but is
% used for scaling of plots.
%
% [Pxx,F] = PSD(X,NFFT,Fs,WINDOW,NOVERLAP) returns a vector of frequen-
% cies the same size as Pxx at which the PSD is estimated, and overlaps
% the sections of X by NOVERLAP samples.
%
% [Pxx, Pxxc, F] = PSD(X,NFFT,Fs,WINDOW,NOVERLAP,P) where P is a scalar
% between 0 and 1, returns the P*100% confidence interval for Pxx.
%
% PSD(X,...,DFLAG), where DFLAG can be 'linear', 'mean' or 'none',
% specifies a detrending mode for the prewindowed sections of X.
% DFLAG can take the place of any parameter in the parameter list
% (besides X) as long as it is last, e.g. PSD(X,'mean');
%
% PSD with no output arguments plots the PSD in the current figure window,
% with confidence intervals if you provide the P parameter.
%
% The default values for the parameters are NFFT = 256 (or LENGTH(X),
% whichever is smaller), NOVERLAP = 0, WINDOW = HANNING(NFFT), Fs = 2,
% P = .95, and DFLAG = 'none'. You can obtain a default parameter by
% leaving it off or inserting an empty matrix [], e.g. PSD(X,[],10000).
%
if nargin < 2, nfft = 512; end
dflag = 'none';
if nargin < 6, LOG = 0; end
if nargin < 5, noverlap = 0; end
if nargin < 4, window = hanning(nfft/2); end
if nargin < 3, Fs = input ('Please enter sampling frequency '); end
if nargin < 2, nfft = length(x); end
% compute PSD
window = window(:);
n = length(x); % Number of data points
nwind = length(window); % length of window
if n < nwind % zero-pad x if it has length less than the window length
x(nwind)=0; n=nwind;
end
x = x(:); % Make sure x is a column vector
% this line must be AFTER the zero-padding in
% case x is a scalar
k = fix((n-noverlap)/(nwind-noverlap)); % Number of windows
% (k = fix(n/nwind) for noverlap=0)
index = 1:nwind;
KMU = k*norm(window)^2; % Normalizing scale factor ==> asymptotically unbiased
Spec = zeros(nfft,1);
for i=1:k
if strcmp(dflag,'none')
xw = window.*(x(index));
elseif strcmp(dflag,'linear')
xw = window.*detrend(x(index));
else
xw = window.*detrend(x(index),0);
end
index = index + (nwind - noverlap);
Xx = abs(fft(xw,nfft)).^2;
Spec = Spec + Xx;
end
% Select first half
if ~any(any(imag(x)~=0)), % if x is not complex
if rem(nfft,2), % nfft odd
select = (1:(nfft+1)/2)';
else
select = (1:nfft/2+1)';
end
Spec = Spec(select)*2;
Spec(1) = Spec(1)/2;
Spec(length(Spec)) = Spec(length(Spec))/2;
else
select = (1:nfft)';
end
freq_vector = (select - 1)*Fs/nfft;
Spec = Spec*(1/KMU); % normalize
f = freq_vector;
if LOG
Spec = 20*log10(Spec);
end