mtmlin.m

多频道功率谱计算程序

function [S,c,Sk]=mtmlin(x,E,V,conf,nfft)
% Syntax: [S,c,Sk]=mtmlin(x,E,V); [S,c,Sk]=mtmlin(x,E,V,conf,nfft)
% Mtmlin produces a linearly weighted Thomson multiple-taper 
% spectral estimate of the time series x, using the
% dpss (tapers) in E and their associated eigenvalues in V.
% Confidence limits (using a chi-squared approach) are
% computed for the level 'conf', which defaults to .95.
% The averaged spectral estimate is returned in S, and the
% confidence limits are in the 2 by N matrix c. The individual
% spectral estimates are returned in Sk. 
% For details, see e. g. Thomson 1982, Park et al. 1987.,
% Percival and Walden 1993.
%
% Written by Eric Breitenberger, version date 10/1/95.
% Please send comments and suggestions to [e-mail]
%

x=x(:);
N=length(x);
W=length(V)/2;
k=2*W-1; % By convention, the first 2W eigenvalues/vectors are stored 

if nargin==3, conf=.95; nfft=2.^(fix(log(N-1)/log(2))+1); end
if nargin==4
  if rem(conf,1), nfft=conf; conf=.95;
  else,  nfft=2.^(fix(log(N-1)/log(2))+1); end
end

% Compute the windowed dfts and the
%   corresponding spectral estimates:
Y=zeros(nfft,k); 
Sk=zeros(nfft,k); 
for i=1:k
  wk=E(:,i).*x;
  Y(:,i)=fft(wk,nfft);
  Sk(:,i)=abs(Y(:,i)).^2;
end

% Select the proper points from fft:
if rem(nfft,2)==0, M=nfft/2+1; else M=(nfft+1)/2; end
%Y=Y(1:M,:);
Sk=Sk(1:M,:);


% Compute the averaged estimate: simple arithmetic
% averaging is used. The Sk can also be weighted 
% by the eigenvalues, as in Park et al. Eqn. 9.;
% note that that eqn. apparently has a typo. as
% the weights should be V and not 1/V.
S=zeros(M,1);
for i=1:k
  S=S+Sk(:,i);
%  S=S+Sk(:,i)*V(i); % Park estimate
end
S=S/k;

% Calculate confidence limits
nu=2*k;
lim=(1-conf)/2;
lim=[lim 1-lim];
conf=wilhil(nu, lim); 
conf=nu./conf;
c(:,1)=conf(1)*S;
c(:,2)=conf(2)*S;