trislice.m

多频道功率谱计算程序

function p=trislice(a,b,lam)
% TRISLICE - Spectrum slicing for a tridiagonal matrix
% Syntax: p=trislice(a,b,lam);
%
% Input:  a - vector of length N, containing the diagonal elements of A.
%         b - vector of length N-1, containing the off-diagonal elements of A.
%       lam - eigenvalue estimate(s). lam may be a vector.
%
% Output: p - number of eigenvalues of A greater than lam. If lam
%              is a vector, pi has n elements, which correspond to the 
%              elements of lam. 
%

a=a(:);
lam=lam(:);
N=length(a);
n=length(lam);
p1=zeros(n,N);
won=ones(n,1);
b=b(:).^2;

p1(:,1)=won*a(1)-lam;

z=find(p1(:,1)==0);
%p1(z,1)=ones(length(z),1)*eps; % replace any zeros with eps.

for i=2:N
  p1(:,i)=won*a(i)-lam-won*b(i-1)./p1(:,i-1);
%  z=find(p1(:,i)==0);
%  p1(z,i)=ones(length(z),i)*eps; % replace any zeros with eps.
end

p=zeros(1,n);
for i=1:n
  p(i)=length(find(p1(i,:)>0));
end