durbin.m

来自「这是个用Durbin递归求预测误差滤波器和他的输出」· M 代码 · 共 36 行

function p = durbin(data,n,q)
%this function is used to sove the toeplitz equation using levinson-durbin
%recrusion. data is the observation time series and n is the length of
%predict-error filter.q is the prediction length. 
%y is the output of the prediction-error filter and a
%is the prediction-error filter.n must be selected so that the length of
%wavelet <= n <length(y).the auto-correlation toeplitz matrix is estimated
%by using FFT techology.
if nargin ~= 3
    error('the input is wrong and we should check it!');
else
    if (n>length(data)||n<=30)
        error('please reselect the n!');
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
l = length(data);
Y = pra_fft(data);
RY = ifft(Y.*conj(Y));
R = (1/n)*real(RY(1:n));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
e = zeros(1,l);
a = zeros(1,l);
e(1) = R(1); 
a(1) = R(2)/R(1);
for k = 2:1:n
    kl = k-1;
    b = sum(R(k:-1:2)'.*a(1:1:kl));
    a(k) = -(R(k)-b)/e(kl);
    a(1:1:kl) = a(1:1:kl) + a(kl:-1:1)*a(k);
    if k~=n
        a = [a(1:kl),a(k)];
       e(k) = e(k-1)*(1-a(k)^2);
   end
end
p = [1,zeros(1,q-1),-a];