Levinson-Durbin算法的实现

    在实现维纳滤波器和预测器的时候，需要计算数据的自相关矩阵的逆。但是当数据量比较大的时候，计算矩阵的逆花费的代价比较大，所以需要使用Levinson-Durbin算法来实现系数的求解。

一、数据模型

k阶前向维纳预测器：

                                    

对上述模型进行变换：

                                    

所以a(0) = 1 , a(i) = -w(i)

注意，在Levinson-Durbin算法中，求解的是a(i), 不是w(i)。

二、Levinson-Durbin迭代算法的实现步骤

    其实关于Levinson-Durbin迭代算法的原理，如果不太清楚的话，也没有关系。只要根据下列步骤就能实现迭代算法。m是预测器的阶数。

1、

2、

3、

4、

初始化值：p(0) = r(0),  a(m-1, 0) = 1,    a(m-1, m) = 0。

其中 r( i ) 是数据的自相关矩阵。

三、matlab实现代码

%alex
clear all;clc
N = 1024;
u = randn(1,N); %产生均值为0，方差为1的高斯白噪声

y(1:2) = u(1:2);
y(3) = -0.69*y(1)+u(3);
y(4) = -0.69*y(2)-0.044*y(1)+u(4);
for n=5:N
    y(n) = -0.69*y(n-2)-0.044*y(n-3)+0.0672*y(n-4)+u(n);
end

%根据耶鲁沃克方程实现如下函数
%求得是前向预测器，所以只用到了数据的自相关。若是滤波则需要用到互相关
%M是预测器的阶数，改变M即可改变预测器的阶数
%ry是训练数据的输入自相关

M = 4;   %M为预测器的阶数
[ ry ] = myxcorr(y,M+1);
Ry = toeplitz(ry(1:M));
Rxy = ry(2:M+1)';
w = inv(Ry)*Rxy;
w = w'

%根据levinson迭代公式实现如下函数
%M是预测器的阶数，改变M即可改变预测器的阶数
%ry是训练数据的输入自相关

%init
p(1) = ry(1);
for m = 1:M
    a(m,1) = 1;
    a(m,m+1) = 0;
    tmp = 0;
    for L = 0:m-1
        tmp = tmp + ry(m-L+1)*a(m-1+1,L+1);   
    end
    delt(m) = tmp;            %delt(1)相当于delta（0）
    k(m) = -delt(m)/p(m);     %k(1)就是看k(1)
    p(m+1) = p(m)*(1-abs(k(m))^2); %p(m+1)相当于p(m)就是p(2)相当于p(1)
    for L = 0:m
        a(m+1,m-L+1) = a(m-1+1,m-L+1) + k(m)*a(m-1+1,L+1);
    end    
end
levinson_w = -a(M+1,2:M+1)
levinson_a = a(M+1,:)

%采用matlab库函数所求的系数
[r,lg] = xcorr(y,'biased');
r(lg<0) = [];
stand_ar = levinson(r,4)


输出结果：

w =
    0.0327   -0.6450   -0.0261    0.0731

levinson_w =
    0.0327   -0.6450   -0.0261    0.0731

levinson_a =
    1.0000   -0.0327    0.6450    0.0261   -0.0731

stand_ar =
    1.0000   -0.0327    0.6441    0.0260   -0.0725

function [ rx ] = myxcorr( x, M )
%该函数用于求出 x 的自相关函数,从r(1)到r(M)
% 输入x为一个1 X N 维矩阵的矩阵，M阶数
% 输出rx为自相关函数
N = length(x);
rx = zeros(1,M);
for i = 1:M    %加一是为了求Rxy
    rx(i) = x(i:N)*x(1:N-i+1)'/(N-i+1);
end
% Rx = toeplitz(rx);
end