pds.m

英文书《Digital Signal Processing with Examples in MATLAB》附带的MATLAB实例

function [P,nsgmts,v]=pds(x,y,N,windo,overlap)
% [P,nsgmts,v]=pds(x,y,N,windo,overlap)
%
% Power density spectrum, averaged over segments of x.
% 
% Inputs:
%  x,y     = Input signal vectors (x,x for auto-spectrum;
%            different vectors for cross-spectrum).
%  N       = segment length (even).  # spectral components=N/2+1.
%  windo   = Data window type:
%                1) Boxcar            4) Hanning
%                2) Tapered           5) Hamming
%                3) Triangular        6) Blackman
%  overlap = Fraction that each data segment of size N
%            overlaps its predecessor.  Must be greater than
%                or equal 0 and less than 1.
% Outputs:
%  P       = Power density spectrum, P(1:N), such that
%            P(n)=P(N-n+2) for n=2,...,N.  If y=x, P =avg.
%            periodogram of x; hence sum(P)/N ~= mean(x.^2).
%  nsgmts  = Number of overlapping segments averaged together.
%  v       = Vector of frequencies. P(n) is at v(n) Hz-s.
% See also power_gain, gain, gain_f

x=row_vec(x);
y=row_vec(y);
L=min(length(x),length(y));

if(N<8 | mod(N,2)==1),
   error('pds: DFT size must be even and at least 8.');
elseif L<N,
   error('pds: Length of x and/or y is < DFT size.');
end

P=zeros(1,N);
nshift=min(N,max(1,round(N*(1-overlap))));
nsgmts=fix(1+(L-N)/nshift);
w=window(N,windo);
for isegmt=0:nsgmts-1
   xx=w.*x((nshift*isegmt+1):(nshift*isegmt+N));
   yy=w.*y((nshift*isegmt+1):(nshift*isegmt+N));
   P=P+conj(fft(xx)).*fft(yy);
end
% Note: division by w*w' includes division by N in 
%       accordance with text eq. 7.24.
P=P/(nsgmts*w*w');
v=[0:N-1]/N;
return