www.pudn.com > yuyinxinhaochulisuanfa.rar > bt.m

function [Pxx,f]=bt(rm,M,N,varargin) 
% Blackman-Tukey algorithm to return the frequency response in Pxx 
%       and frequency vector in f. 
% 
%  [Pxx,f]=BT(rm,M,N,'half')  or 
%  [Pxx,f]=BT(rm,M,N,'whole')  or 
%  [Pxx,f]=BT(rm,M,N)  
%         Pxx <--- N-point complex frequency response 
%         f   <--- N-point frequency vector in radians/sample  
%         rm  ---> autocorrelation data from negative lag to positive lag 
%         M   ---> used data of the autocorrelation from -M lag to M lag, 
%                  which should be less than the length of rm 
%         N   ---> the number of computered points in [0,2*pi) 
%     'whole' ---> return the all points in [0,2*pi) 
%      'half' ---> return the half points in [0,pi), which is default . 
% 
%  The formulation for the algorithm is: 
%          M 
%   P(w)= sum  r(m)*exp(-j*w*m),   |M|<=N-1 
%         m=-M       
 
%  Copyright by the Author: Zhilin Zhang 
%  $ Revision 1.5 $$ Date: 2003/04/02 01:39 $ 
% 
%  References 
%      胡广书，数字信号处理——理论、算法与实现》P325，清华大学出版社，1996 
 
 
% test input arguments number 
if nargin<3  
    error('Too few input arguments'); 
elseif nargin>4 
    error('Too many input arguments'); 
end 
 
rlen=length(rm);                     % length of input correlation data 
cent=(rlen+1)/2;                     % locate the center point in the data 
 
% cheke wether r(m) is from negative lag to positive lag 
if rm(cent)<=rm(1) 
    error('The r(m) should be from negative lag to positive lag'); 
end 
 
%%%%%%%%%%%%%%%%%%%%% Initialization %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
fdelt=2*pi/N;                      %  frequency value between neibor points 
 
%  Pxx is N-point frequency response when the fourth arguments is 'whole'; 
Pxx=zeros(1,N);  
 
%%%%%%%%%%%%%%%%%%%% computer %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
for index=0:N-1 
    temp=0; 
    for m=-M:M 
        temp=temp+rm(m+cent)*exp(-j*m*index*fdelt); 
    end 
    Pxx(index+1)=temp; 
end 
f=[0:N-1]*fdelt; 
 
%%%%%%%%%%%%%%%%%%% determine output %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% if the fourth argument is 1, or the fourth argument is NOT determined, 
% then truncation should be apply to the output arguments f and Pxx 
f2=f(1:floor((N+1)/2));                           % truncation of f         
Pxx2=Pxx(1:floor((N+1)/2));                       % truncation of Pxx 
 
if nargin==3 | strcmp(varargin{1},'half') 
    Pxx=Pxx2; 
    f=f2; 
end