burg.m

来自「在MATLAB上实现burg算法」· M 代码 · 共 47 行

function a = burg(x,ORDER)
% MADE BY H.X.GUAN,2008/12/18
% BURG ALGORITHEM,AR model

%%%%% Initial %%%%%%
N = length(x);
f = x; % 前向
g = x; % 后向
% 假设x为行向量
P0 = x*x'/N;
PM = P0;
m = 2;
if ORDER>N
    disp('阶数超过向量长度无法计算');
    a = 0;
    return;
end
a = zeros(1,ORDER+1);
am = a;
a(1) = 1;
am(1) = 1;
%%%%%%%%%%%%%%%%%%%%
while ( m <= ORDER+1 )
    sum1 = 0;
    sum2 = 0;
    for n=m:N
        sum1 = sum1 + f(n)*g(n-1);
        sum2 = sum2 + f(n)*f(n) + g(n-1)*g(n-1);
    end
    Km = -1.0*sum1 / (sum2*0.5);
    for i=2:m
        am(i) = a(i) + Km*a(m-i+1);
    end
%    am(m) = Km;
    for i=2:m
        a(i) = am(i);
    end
    PM = (1-Km*Km)*PM;
    for n=2:N
        fm(n) = f(n) + Km*g(n-1);
        gm(n) = Km*f(n) + g(n-1);
    end
    f = fm;
    g = gm;
    m = m + 1;
end