dpssint.m

多频道功率谱计算程序

function [En,V] = dpssint(N, W, M, int) 
% Syntax: [En,V]=dpssint(N,W); [En,V]=dpssint(N,W,M,'spline'); 
%  Dpssint calculates discrete prolate spheroidal 
%  sequences for the parameters N and W. Note that 
%  W is normally 2, 5/2, 3, 7/2, or 4 - not i/N. The  
%  dpss are interpolated from previously calculated  
%  dpss of order M (128, 256, 512, or 1024). 256 is the  
%  default for M. The interpolation can be 'linear'  
%  or 'spline'. 'Linear' is faster, 'spline' the default. 
%  Linear interpolation can only be used for M>N.  
%  Returns: 
%              E: matrix of dpss (N by 2W) 
%              V: eigenvalue vector (2W) 
%  
% Errors in the interpolated dpss are very small but should be  
% checked if possible. The differences between interpolated 
% values and values from dpsscalc are generally of order 
% 10ee-5 or better. Spline interpolation is generally more 
% accurate. Fractional errors can be very large near 
% the zero-crossings but this does not seriously affect 
% windowing calculations. The error can be reduced by using 
% values for M which are close to N. 
% 
% Written by Eric Breitenberger, version date 10/3/95. 
% Please send comments and suggestions to [e-mail] 
% 
 
if     nargin==2, 
  M=256; int='spline'; 
elseif nargin==3, 
  if isstr(M), int=M; M=256;  
  else, int='spline';, end 
end 
 
if int=='linear' & N>M 
  error('Linear interpolation cannot be used for N>M. Use splining instead.') 
end 
 
if rem(W,1)==0 
  filename=['s' int2str(M) '_' int2str(W) '.mat'];  
else 
  filename=['s' int2str(M) '_' int2str(10*W) '.mat'];  
end 
 
if exist(filename)==2 
  eval(['load ' filename]) 
else 
  filename 
  error('Requested DPSS file not found.') 
end 
 
 
k=2*W;  % Use only first k values 
E=E(:,1:k); 
V=V(1:k); 
En=zeros(N,k); 
x=1:M; 
 
% The scaling for the interpolation: 
% This is not necessarily optimal, and  
% changing s can improve accuracy. 
  
s=M/N; 
midm=(M+1)/2; 
midn=(N+1)/2; 
delta=midm-s*midn; 
xi=linspace(1-delta, M+delta, N); 
 
% Interpolate from M values to N 
% Spline interpolation is a bit better, 
% but takes about twice as long. 
% Errors from linear interpolation are  
% usually smaller than errors from scaling. 
 
if int=='linear' 
  for i=1:k, En(:,i)=interp1(x,E(:,i),xi); end 
elseif int=='spline' 
  for i=1:k, En(:,i)=spline(x,E(:,i),xi)'; end 
end 
 
% Re-normalize the eigenvectors 
En=En./(ones(N,1)*sqrt(sum(En.*En)));