二阶格型滤波器matlab的实现.m

格型滤波器(lattice filter)Matlab的实现帖子,只用于2阶。请高手看看。不过我还是不会编更高阶的格型滤波器。欢迎一起探讨！[e-mail]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   Lattice_AF.m
%   Lattice adaptive filter
%   LMS算法格型滤波器Matlab的实现帖子！！只用于2阶。请高手看看！！
%   不过我还是不会编更高阶的格型滤波器！！欢迎一起探讨！[e-mail]
%   2006.06.30
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
close all
hold off
%channel system order
sysorder = 2 ;
snr=30 ;
% Number of system points
N=2000;
inp = randn(N,1);
n = randn(N,1)/10.^(snr/10);
b=[1,0,0]; 
a=[1,-1.558,0.81];
Gz = tf(b,a,1);
h= [1.558;-0.81];
y = lsim(Gz,inp);% inp通过gz的输出响应
%add some noise
d = y + n;
totallength=size(d,1);
%begin of algorithm
M=totallength;
%初始化各系数
f0=zeros(1,M);
f1=zeros(1,M);
f2=zeros(1,M);
k1=zeros(1,M);
k2=zeros(1,M);
b0=zeros(1,M);
b1=zeros(1,M);
b2=zeros(1,M);
a1=zeros(1,M);
a2=zeros(1,M);
f0=inp;
b0=inp;
f1(1)=inp(1);
k1(1)=0;
k2(1)=0;
b1(1)=0;
k1(2)=0;
f2(1)=inp(1);
b2(1)=0;
k2(2)=0;
u=0.05;
for n = sysorder : M-1
    i=n+2-sysorder;
    f1(i)=f0(i)+k1(i)*b0(i-1);
    b1(i)=b0(i-1)+k1(i)*f0(i);
    k1(i+1)=k1(i)-2*u*(f1(i)*b0(i-1)+b1(i)*f0(i));
    f2(i)=f1(i)+k2(i)*b1(i-1);
    b2(i)=b1(i-1)+k2(i)*f1(i);
    k2(i+1)=k2(i)-2*u*(f2(i)*b1(i-1)+b2(i)*f1(i));
end 
for n = sysorder : M
    i=n+1-sysorder;
    a1(i)=(-k1(i)*(1+k2(i)));
    a2(i)=-k2(i);
end
hold on
figure(1);
subplot(3,1,1)
plot(d)
plot(y,'r');
title('System output') ;
xlabel('Samples')
ylabel('True and estimated output')
subplot(3,1,3)
plot(h, 'k+')
hold on
plot(a, 'r*')
legend('Actual weights','Estimated weights')
title('Comparison of the actual weights and the estimated weights') ;
axis([0 3 -2.35 2.35])