ld.m

采用matlab编写的Levinson Durbin 算法

%Levinson Durbin algorithm to spectra, reference function "pyulear"
%pruduce the signal

close all;
clear all;
N = 1000;% Number of system points
input = wgn(N,1,0);
A = [0, 0.8];%b0 = 1,a0 = 1,a1 = 0, a2 = 0.8
x = filter(1,[1,A],input);
% x = filter([1,1,1],1,input);
%zplane(1,[1,A);
%figure;

%computer the self-R matrix
R = xcorr(x,'coeff') ;
R = R(N:end);
 
p = 2;
a = zeros(p+1,p+1);
sum1 = 0;
sum2 = 0;

for k = 0:p - 1 %computer form k = 0 ?
    a(k + 1, 1) = 1;
    sum1 = 0;
    sum2 = 0;
    for i = 0:k 
        sum1 = sum1 + a(k+1 , i+1) * R(i+1);
    end
    for i = 0:k   
        sum2 = sum2 + a(k+1, i+1) * R(k + 2 - i);        
    end
   % sita_sqr(k) = sum1 + R(1);
    sita_sqr(k+1) = sum1;
    D(k+1) = sum2;
    gamma( k + 2) = D(k+1)/sita_sqr(k+1);
    sita_sqr( k + 2) = (1 - gamma(k + 2)^2) * sita_sqr(k+1);
    for i = 1:k 
         a( k + 2 ,i+1) = a( k+1 , i+1) - gamma( k + 2) * a( k+1 , k + 2 - i);
    end
    a(k + 2,k + 2) = - gamma( k + 2);
end

% for k = 1:p %computer form k = 0 ?
%     a(k, 1) = 1;
%     sum1 = 0;
%     sum2 = 0;
%     for i = 1:k 
%         sum1 = sum1 + a(k , i) * R(i);
%     end
%     for i = 1:k   
%         sum2 = sum2 + a(k, i) * R(k + 2- i);        
%     end
%    % sita_sqr(k) = sum1 + R(1);
%    sita_sqr(k) = sum1;
%     D(k) = sum2;
%     gamma( k + 1) = D(k)/sita_sqr(k);
%     sita_sqr( k + 1) = (1 - gamma(k + 1)^2) * sita_sqr(k);
%     for i = 2:k 
%          a( k + 1 ,i) = a( k , i) - gamma( k + 1) * a( k , k + 2 - i);
%     end
%     a(k + 1,k + 1) = - gamma( k + 1);
% end

M = 200;

for k = 1:M
    sum = 0;
    for i = 2:p+1
        sum = sum + a(p+1,i) * exp( -j * 2 * pi/M*(i-1)*(k-1)) ;
    end
    S(k) = sita_sqr(p+1)/(  abs( 1 + sum )^2 );
    
end
Q = log (S);
plot(S);
figure;
plot(Q);