ak.m

来自「格型滤波器的K值估计以及其与相应线性滤波器之间系数A的转换」· M 代码 · 共 28 行

function ak=AK(u,N)
v=randn(1,N);
x(1)=v(1);
x(2)=1.558*x(1)+v(2);
k1(1)=0;
k2(1)=0;
k1(2)=0;
k2(2)=0;
for n=3:N
    x(n)=1.558*x(n-1)-0.81*x(n-2)+v(n);
    e0f(n-1)=x(n-1);
    e0b(n-2)=x(n-2);
    e1f(n-1)=x(n-1)+k1(n-1)*x(n-2);
    e1b(n-1)=x(n-2)+k1(n-1)*x(n-1);
    if n==3
        e1b(n-2)=k1(n-2)*x(n-2);
    else
    e1b(n-2)=x(n-3)+k1(n-2)*x(n-2);
    end
    e2f(n-1)=e1f(n-1)+k2(n-1)*e1b(n-2);
    e2b(n-1)=e1b(n-2)+k2(n-1)*e1f(n-1);
    k1(n)=k1(n-1)-2*u*[e1f(n-1)*e0b(n-2)+e1b(n-1)*e0f(n-1)];
    k2(n)=k2(n-1)-2*u*[e2f(n-1)*e1b(n-2)+e2b(n-1)*e1f(n-1)];
    a1(n)=-k1(n)*(1+k2(n));
    a2(n)=-k2(n);
end
t=1:N;
plot(t,a1,'r',t,a2)