427 lines
14 KiB
Matlab
427 lines
14 KiB
Matlab
function [Spec,f,SpecConf] = spectrum2(varargin)
|
|
%SPECTRUM Power spectrum estimate of one or two data sequences.
|
|
% P=SPECTRUM2(X,NFFT,NOVERLAP,WIND) estimates the Power Spectral Density of
|
|
% signal vector(s) X using Welch's averaged periodogram method. X must be a
|
|
% set of column vectors. They are divided into overlapping sections, each of
|
|
% which is detrended and windowed by the WINDOW parameter, then zero padded
|
|
% to length NFFT. The magnitude squared of the length NFFT DFTs of the sec-
|
|
% tions are averaged to form Pxx. P is a two times SIZE(X,2) column matrix
|
|
% P = [Pxx Pxxc]; the second half of columns Pxxc are the 95% confidence
|
|
% intervals. The number of rows of P is 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.
|
|
%
|
|
% [P,F] = SPECTRUM2(X,NFFT,NOVERLAP,WINDOW,Fs) given a sampling frequency
|
|
% Fs returns a vector of frequencies the same length as Pxx at which the
|
|
% PSD is estimated. PLOT(F,P(:,1:end/2)) plots the power spectrum estimate
|
|
% versus true frequency.
|
|
%
|
|
% [P, F] = SPECTRUM2(X,NFFT,NOVERLAP,WINDOW,Fs,Pr) where Pr is a scalar
|
|
% between 0 and 1, overrides the default 95% confidence interval and
|
|
% returns the Pr*100% confidence interval for Pxx instead.
|
|
%
|
|
% SPECTRUM2(X) with no output arguments plots the PSD in the current
|
|
% figure window, with confidence intervals.
|
|
%
|
|
% The default values for the parameters are NFFT = 256 (or SIZE(X,1),
|
|
% whichever is smaller), NOVERLAP = 0, WINDOW = HANNING(NFFT), Fs = 2,
|
|
% and Pr = .95. You can obtain a default parameter by leaving it out
|
|
% or inserting an empty matrix [], e.g. SPECTRUM(X,[],128).
|
|
%
|
|
% P = SPECTRUM2(X,Y) performs spectral analysis of the two times SIZE(X,2)
|
|
% sequences X and Y using the Welch method. SPECTRUM returns the 8 times
|
|
% SIZE(X,2)column array
|
|
% P = [Pxx Pyy Pxy Txy Cxy Pxxc Pyyc Pxyc]
|
|
% where
|
|
% Pxx = X-vector power spectral density
|
|
% Pyy = Y-vector power spectral density
|
|
% Pxy = Cross spectral density
|
|
% Txy = Complex transfer function from X to Y = Pxy./Pxx
|
|
% Cxy = Coherence function between X and Y = (abs(Pxy).^2)./(Pxx.*Pyy)
|
|
% Pxxc,Pyyc,Pxyc = Confidence range.
|
|
% All input and output options are otherwise exactly the same as for the
|
|
% single input case.
|
|
%
|
|
% SPECTRUM2(X,Y) with no output arguments will plot Pxx, Pyy, abs(Txy),
|
|
% angle(Txy) and Cxy in sequence, pausing between plots.
|
|
%
|
|
% SPECTRUM2(X,...,DFLAG), where DFLAG can be 'linear', 'mean' or 'none',
|
|
% specifies a detrending mode for the prewindowed sections of X (and Y).
|
|
% DFLAG can take the place of any parameter in the parameter list
|
|
% (besides X) as long as it is last, e.g. SPECTRUM(X,'none');
|
|
%
|
|
% See also SPECTRUM, PSD, CSD, TFE, COHERE, SPECGRAM, SPECPLOT, DETREND,
|
|
% PMTM, PMUSIC.
|
|
% ETFE, SPA, and ARX in the Identification Toolbox.
|
|
|
|
% The units on the power spectra Pxx and Pyy are such that, using
|
|
% Parseval's theorem:
|
|
%
|
|
% SUM(Pxx)/LENGTH(Pxx) = SUM(X.^2)/size(x,1) = COV(X)
|
|
%
|
|
% The RMS value of the signal is the square root of this.
|
|
% If the input signal is in Volts as a function of time, then
|
|
% the units on Pxx are Volts^2*seconds = Volt^2/Hz.
|
|
%
|
|
% Here are the covariance, RMS, and spectral amplitude values of
|
|
% some common functions:
|
|
% Function Cov=SUM(Pxx)/LENGTH(Pxx) RMS Pxx
|
|
% a*sin(w*t) a^2/2 a/sqrt(2) a^2*LENGTH(Pxx)/4
|
|
%Normal: a*rand(t) a^2 a a^2
|
|
%Uniform: a*rand(t) a^2/12 a/sqrt(12) a^2/12
|
|
%
|
|
% For example, a pure sine wave with amplitude A has an RMS value
|
|
% of A/sqrt(2), so A = SQRT(2*SUM(Pxx)/LENGTH(Pxx)).
|
|
%
|
|
% See Page 556, A.V. Oppenheim and R.W. Schafer, Digital Signal
|
|
% Processing, Prentice-Hall, 1975.
|
|
|
|
error(nargchk(1,8,nargin))
|
|
[msg,x,y,nfft,noverlap,window,Fs,p,dflag]=specchk2(varargin);
|
|
error(msg)
|
|
|
|
|
|
if isempty(p),
|
|
p = .95; % default confidence interval even if not asked for
|
|
end
|
|
|
|
n = size(x,1); % Number of data points
|
|
ns = size(x,2); % Number of signals
|
|
nwind = length(window);
|
|
if n < nwind % zero-pad x (and y) if length less than the window length
|
|
x(nwind,:)=0; n=nwind;
|
|
if ~isempty(y), y(nwind)=0; end
|
|
end
|
|
|
|
k = fix((n-noverlap)/(nwind-noverlap)); % Number of windows
|
|
index = 1:nwind;
|
|
KMU = k*norm(window)^2; % Normalizing scale factor ==> asymptotically unbiased
|
|
% KMU = k*sum(window)^2;% alt. Nrmlzng scale factor ==> peaks are about right
|
|
|
|
if (isempty(y)) % Single sequence case.
|
|
Pxx = zeros(nfft,ns); Pxx2 = zeros(nfft,ns);
|
|
for i=1:k
|
|
for l=1:ns
|
|
if strcmp(dflag,'linear')
|
|
xw = window.*detrend(x(index,l));
|
|
elseif strcmp(dflag,'none')
|
|
xw = window.*(x(index,l));
|
|
else
|
|
xw = window.*detrend(x(index,l),0);
|
|
end
|
|
Xx = abs(fft(xw,nfft)).^2;
|
|
Pxx(:,l) = Pxx(:,l) + Xx;
|
|
Pxx2(:,l) = Pxx2(:,l) + abs(Xx).^2;
|
|
end
|
|
index = index + (nwind - noverlap);
|
|
end
|
|
% Select first half
|
|
if ~any(any(imag(x)~=0)), % if x and y are not complex
|
|
if rem(nfft,2), % nfft odd
|
|
select = [1:(nfft+1)/2];
|
|
else
|
|
select = [1:nfft/2+1]; % include DC AND Nyquist
|
|
end
|
|
else
|
|
select = 1:nfft;
|
|
end
|
|
Pxx = Pxx(select,:);
|
|
Pxx2 = Pxx2(select,:);
|
|
cPxx = zeros(size(Pxx));
|
|
if k > 1
|
|
c = (k.*Pxx2-abs(Pxx).^2)./(k-1);
|
|
c = max(c,zeros(size(Pxx)));
|
|
cPxx = sqrt(c);
|
|
end
|
|
ff = sqrt(2)*erfinv(p); % Equal-tails.
|
|
Pxxc = ff.*cPxx/KMU;
|
|
P = Pxx/KMU;
|
|
Pc = Pxxc;
|
|
else
|
|
Pxx = zeros(nfft,ns); % Dual sequence case.
|
|
Pxy = Pxx; Pxx2 = Pxx; Pxy2 = Pxx;
|
|
Pyy = zeros(nfft,1); Pyy2 = Pyy;
|
|
for i=1:k
|
|
if strcmp(dflag,'linear')
|
|
yw = window.*detrend(y(index));
|
|
elseif strcmp(dflag,'none')
|
|
yw = window.*(y(index));
|
|
else
|
|
yw = window.*detrend(y(index),0);
|
|
end
|
|
Yy = fft(yw,nfft);
|
|
Yy2 = abs(Yy).^2;
|
|
Pyy = Pyy + Yy2;
|
|
Pyy2 = Pyy2 + abs(Yy2).^2;
|
|
for l=1:ns
|
|
if strcmp(dflag,'linear')
|
|
xw = window.*detrend(x(index,l));
|
|
elseif strcmp(dflag,'none')
|
|
xw = window.*(x(index,l));
|
|
else
|
|
xw = window.*detrend(x(index,l),0);
|
|
end
|
|
Xx = fft(xw,nfft);
|
|
Xx2 = abs(Xx).^2;
|
|
Pxx(:,l) = Pxx(:,l) + Xx2;
|
|
Pxx2(:,l) = Pxx2(:,l) + abs(Xx2).^2;
|
|
Xy = Yy .* conj(Xx);
|
|
Pxy(:,l) = Pxy(:,l) + Xy;
|
|
Pxy2(:,l) = Pxy2(:,l) + Xy .* conj(Xy);
|
|
end
|
|
index = index + (nwind - noverlap);
|
|
end
|
|
% Select first half
|
|
if ~any(any(imag([x y])~=0)), % if x and y are not complex
|
|
if rem(nfft,2), % nfft odd
|
|
select = [1:(nfft+1)/2];
|
|
else
|
|
select = [1:nfft/2+1]; % include DC AND Nyquist
|
|
end
|
|
else
|
|
select = 1:nfft;
|
|
end
|
|
Pxx = Pxx(select,:);
|
|
Pxy = Pxy(select,:);
|
|
Pxx2 = Pxx2(select,:);
|
|
Pxy2 = Pxy2(select,:);
|
|
Pyy = Pyy(select);
|
|
Pyy2 = Pyy2(select);
|
|
cPxx = zeros(size(Pxx));
|
|
cPyy = zeros(size(Pyy));
|
|
cPxy = cPxx;
|
|
if k > 1
|
|
c = max((k.*Pxx2-abs(Pxx).^2)./(k-1),zeros(size(Pxx)));
|
|
cPxx = sqrt(c);
|
|
c = max((k.*Pyy2-abs(Pyy).^2)./(k-1),zeros(size(Pyy)));
|
|
cPyy = sqrt(c);
|
|
c = max((k.*Pxy2-abs(Pxy).^2)./(k-1),zeros(size(Pxx)));
|
|
cPxy = sqrt(c);
|
|
end
|
|
Txy = Pxy./Pxx;
|
|
Cxy = (abs(Pxy).^2)./(Pxx.*repmat(Pyy,1,size(Pxx,2)));
|
|
ff = sqrt(2)*erfinv(p); % Equal-tails.
|
|
Pxx = Pxx/KMU;
|
|
Pyy = Pyy/KMU;
|
|
Pxy = Pxy/KMU;
|
|
Pxxc = ff.*cPxx/KMU;
|
|
Pxyc = ff.*cPxy/KMU;
|
|
Pyyc = ff.*cPyy/KMU;
|
|
P = [Pxx Pyy Pxy Txy Cxy];
|
|
Pc = [Pxxc Pyyc Pxyc];
|
|
end
|
|
freq_vector = (select - 1)'*Fs/nfft;
|
|
|
|
if nargout == 0, % do plots
|
|
newplot;
|
|
if Fs==2, xl='Frequency'; else, xl = 'f [Hz]'; end
|
|
nplot=1+(1-isempty(y))*4;
|
|
|
|
subplot(nplot,1,1);
|
|
c = [max(Pxx-Pxxc,0) Pxx+Pxxc];
|
|
c = c.*(c>0);
|
|
|
|
h=semilogy(freq_vector,Pxx,...
|
|
freq_vector,c(:,1:size(c,2)/2),'--',...
|
|
freq_vector,c(:,size(c,2)/2+1:end),'--');
|
|
title('\bf X Power Spectral Density')
|
|
ylabel('P_x')
|
|
xlabel(xl)
|
|
if length(h)>3
|
|
s={};
|
|
for k=1:length(h)/3
|
|
c=get(h(k),'Color');
|
|
set(h(k+length(h)/3),'Color',c);
|
|
set(h(k+2*length(h)/3),'Color',c);
|
|
s{k}=['x_' num2str(k)];
|
|
end
|
|
legend(h(1:length(h)/3),s);
|
|
end
|
|
if (isempty(y)), % single sequence case
|
|
return
|
|
end
|
|
|
|
subplot(nplot,1,2);
|
|
c = [max(Pyy-Pyyc,0) Pyy+Pyyc];
|
|
c = c.*(c>0);
|
|
h=semilogy(freq_vector,Pyy,...
|
|
freq_vector,c(:,1),'--',...
|
|
freq_vector,c(:,2),'--');
|
|
if size(Pxx,2)>1
|
|
for k=1:length(h)/3
|
|
c=get(h(k),'Color');
|
|
set(h(k+length(h)/3),'Color',c);
|
|
set(h(k+2*length(h)/3),'Color',c);
|
|
end
|
|
end
|
|
|
|
title('\bf Y Power Spectral Density')
|
|
ylabel('P_y')
|
|
xlabel(xl)
|
|
|
|
subplot(nplot,1,3);
|
|
semilogy(freq_vector,abs(Txy));
|
|
title('\bf Transfer function magnitude')
|
|
ylabel('T_{xy}^{(m)}')
|
|
xlabel(xl)
|
|
|
|
subplot(nplot,1,4);
|
|
plot(freq_vector,(angle(Txy))), ...
|
|
title('\bf Transfer function phase')
|
|
ylabel('T_{xy}^{(p)}')
|
|
xlabel(xl)
|
|
|
|
subplot(nplot,1,5);
|
|
plot(freq_vector,Cxy);
|
|
title('\bf Coherence')
|
|
ylabel('C_{xy}')
|
|
xlabel(xl)
|
|
if exist ('niceaxes') ~= 0
|
|
niceaxes(findall(gcf,'Type','axes'));
|
|
end
|
|
elseif nargout ==1,
|
|
Spec = P;
|
|
elseif nargout ==2,
|
|
Spec = P;
|
|
f = freq_vector;
|
|
elseif nargout ==3,
|
|
Spec = P;
|
|
f = freq_vector;
|
|
SpecConf = Pc;
|
|
end
|
|
|
|
function [msg,x,y,nfft,noverlap,window,Fs,p,dflag] = specchk2(P)
|
|
%SPECCHK Helper function for SPECTRUM
|
|
% SPECCHK(P) takes the cell array P and uses each cell as
|
|
% an input argument. Assumes P has between 1 and 7 elements.
|
|
|
|
msg = [];
|
|
if size(P{1},1)<=1
|
|
if max(size(P{1}))==1
|
|
msg = 'Input data must be a vector, not a scalar.';
|
|
else
|
|
msg = 'Requires column vector input.';
|
|
end
|
|
x=[];
|
|
y=[];
|
|
elseif (length(P)>1),
|
|
if (all(size(P{1},1)==size(P{2})) & (size(P{1},1)>1) ) | ...
|
|
size(P{2},1)>1, % 0ne signal or 2 present?
|
|
% two signals, x and y, present
|
|
x = P{1}; y = P{2};
|
|
% shift parameters one left
|
|
P(1) = [];
|
|
else
|
|
% only one signal, x, present
|
|
x = P{1}; y = [];
|
|
end
|
|
else % length(P) == 1
|
|
% only one signal, x, present
|
|
x = P{1}; y = [];
|
|
end
|
|
|
|
% now x and y are defined; let's get the rest
|
|
|
|
if length(P) == 1
|
|
nfft = min(size(x,1),256);
|
|
window = hanning(nfft);
|
|
noverlap = 0;
|
|
Fs = 2;
|
|
p = [];
|
|
dflag = 'linear';
|
|
elseif length(P) == 2
|
|
if isempty(P{2}), dflag = 'linear'; nfft = min(size(x,1),256);
|
|
elseif isstr(P{2}), dflag = P{2}; nfft = min(size(x,1),256);
|
|
else dflag = 'linear'; nfft = P{2}; end
|
|
window = hanning(nfft);
|
|
noverlap = 0;
|
|
Fs = 2;
|
|
p = [];
|
|
elseif length(P) == 3
|
|
if isempty(P{2}), nfft = min(size(x,1),256); else nfft=P{2}; end
|
|
if isempty(P{3}), dflag = 'linear'; noverlap = 0;
|
|
elseif isstr(P{3}), dflag = P{3}; noverlap = 0;
|
|
else dflag = 'linear'; noverlap = P{3}; end
|
|
window = hanning(nfft);
|
|
Fs = 2;
|
|
p = [];
|
|
elseif length(P) == 4
|
|
if isempty(P{2}), nfft = min(size(x,1),256); else nfft=P{2}; end
|
|
if isstr(P{4})
|
|
dflag = P{4};
|
|
window = hanning(nfft);
|
|
else
|
|
dflag = 'linear';
|
|
window = P{4}; window = window(:); % force window to be a column
|
|
if length(window) == 1, window = hanning(window); end
|
|
if isempty(window), window = hanning(nfft); end
|
|
end
|
|
if isempty(P{3}), noverlap = 0; else noverlap=P{3}; end
|
|
Fs = 2;
|
|
p = [];
|
|
elseif length(P) == 5
|
|
if isempty(P{2}), nfft = min(size(x,1),256); else nfft=P{2}; end
|
|
window = P{4}; window = window(:); % force window to be a column
|
|
if length(window) == 1, window = hanning(window); end
|
|
if isempty(window), window = hanning(nfft); end
|
|
if isempty(P{3}), noverlap = 0; else noverlap=P{3}; end
|
|
if isstr(P{5})
|
|
dflag = P{5};
|
|
Fs = 2;
|
|
else
|
|
dflag = 'linear';
|
|
if isempty(P{5}), Fs = 2; else Fs = P{5}; end
|
|
end
|
|
p = [];
|
|
elseif length(P) == 6
|
|
if isempty(P{2}), nfft = min(size(x,1),256); else nfft=P{2}; end
|
|
window = P{4}; window = window(:); % force window to be a column
|
|
if length(window) == 1, window = hanning(window); end
|
|
if isempty(window), window = hanning(nfft); end
|
|
if isempty(P{3}), noverlap = 0; else noverlap=P{3}; end
|
|
if isempty(P{5}), Fs = 2; else Fs = P{5}; end
|
|
if isstr(P{6})
|
|
dflag = P{6};
|
|
p = [];
|
|
else
|
|
dflag = 'linear';
|
|
if isempty(P{6}), p = .95; else p = P{6}; end
|
|
end
|
|
elseif length(P) == 7
|
|
if isempty(P{2}), nfft = min(size(x,1),256); else nfft=P{2}; end
|
|
window = P{4}; window = window(:); % force window to be a column
|
|
if length(window) == 1, window = hanning(window); end
|
|
if isempty(window), window = hanning(nfft); end
|
|
if isempty(P{3}), noverlap = 0; else noverlap=P{3}; end
|
|
if isempty(P{5}), Fs = 2; else Fs = P{5}; end
|
|
if isempty(P{6}), p = .95; else p = P{6}; end
|
|
if isstr(P{7})
|
|
dflag = P{7};
|
|
else
|
|
msg = 'DFLAG parameter must be a string.'; return
|
|
end
|
|
end
|
|
|
|
% NOW do error checking
|
|
if isempty(msg)
|
|
if (nfft<length(window)),
|
|
msg = 'Requires window''s length to be no greater than the FFT length.';
|
|
end
|
|
if (noverlap >= length(window)),
|
|
msg = 'Requires NOVERLAP to be strictly less than the window length.';
|
|
end
|
|
if (nfft ~= abs(round(nfft)))|(noverlap ~= abs(round(noverlap))),
|
|
msg = 'Requires positive integer values for NFFT and NOVERLAP.';
|
|
end
|
|
if ~isempty(p),
|
|
if (prod(size(p))>1)|(p(1,1)>1)|(p(1,1)<0),
|
|
msg = 'Requires confidence parameter to be a scalar between 0 and 1.';
|
|
end
|
|
end
|
|
if (min(size(y))~=1)&(~isempty(y)),
|
|
msg = 'Requires column vector input as second signal.';
|
|
end
|
|
if (size(x,1)~=length(y))&(~isempty(y)),
|
|
msg = 'Requires X and Y to have the same number of rows.';
|
|
end
|
|
end |