dpsscal2.m

多频道功率谱计算程序

  function [E,V] = dpsscal2(N,W)
% DPCSSCAL2 - Slepian Sequence calculation for short sequences
%  Syntax:  [E,V] = dpsscal2(N,W);
%
%  Dpsscal2 calculates discrete prolate spheroidal sequences for the
%  parameters N and W. The tridiagonal method of Slepian is used.
%  Note that W is of form 2, 5/2, 3, 7/2, 4, ... and not 2/N, 5/2N, 3/N, etc.
%
%  DPSSCAL2 is very inefficient computationally for large N. If N is greater
%  than about 200 or so (this depends on the platform), DPSSCALC should
%  be used instead. 
%
%  Returns:
%              E: matrix of dpss (N by 2W)
%              V: eigenvalue vector (2W)
%
% See Percival and Walden, Chapter 8.
%
% Written by Eric Breitenberger, version date 2/6/95.
% Please send comments and suggestions to [e-mail]
%

W=W/N;
% Generate the diagonal
d=((N-1-2*(0:N-1)).^2)*.25*cos(2*pi*W);

% Generate the off-diagonals
off=(1:N-1).*(N-1:-1:1)/2;

% Generate the tri-diagonal matrix
B = diag(d) + diag(off,1) + diag(off, -1);

% Get the eigenvectors/values of B, then sort
[E, V]=eig(B);
[v, index]=sort(diag(V));
clear B V v
E=E(:,index(N:-1:1));
% Keep only the needed eigenvectors
E=E(:,1:2*N*W);

% Now calculate the desired eigenvalues by plugging
% the eigenvectors back into the defining equation:
d(1)=2*W;
d(2:N)=sin(2*pi*W*(1:N-1))./(pi*(1:N-1));

A=toeplitz(d);
V=E'*A*E;
V=diag(V);

% Normalize the eigenvectors
E=E./(ones(N,1)*sqrt(sum(E.*E)));

% Polarize symmetric dpss
d=mean(E);
for i=1:2:2*N*W
  if d(i)<0, E(:,i)=-E(:,i); end
end

% Polarize anti-symmetric dpss
for i=2:2:2*N*W
  if E(2,i)<0, E(:,i)=-E(:,i); end
end