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