clear;
N=200;
order=2;
M=order+1;
t=1;                                %t=δ
%n=1時刻的初始化
A=randn(1,N);                       %噪聲                           
x=zeros(1,N);                       %輸入
y=zeros(1,N);                       %輸出
e=zeros(1,N);
error=zeros(1,N);                   %誤差
eb=zeros(M,N);                      %后向預測誤差
ef=zeros(M,N);                      %前向預測誤差
delta=zeros(M,N);                   %△
gama=ones(M,N);                     %角參量
Ef=t*ones(M,N);                     %前向剩余誤差
Eb=t*ones(M,N);                     %后向剩余誤差
kf=zeros(M-1,N);                    %前向反射系數
kb=zeros(M-1,N);                    %后向反射系數
a=zeros(M-1,N);                     %權矢量系數
%求出輸入x(n)
x(1)=A(1);
x(2)=A(2)+1.558*A(1);
for i=3:N
    x(i)=A(i)+1.558*x(i-1)-0.81*x(i-2);    
end
z1=ones(1,N)*1.58;
z2=ones(1,N)*(-0.81);
%m=1的初始化
eb(1,:)=x;
ef(1,:)=x;
for j=2:N
    Eb(1,j)=Ef(1,j-1)+x(j)*x(j);
    Ef(1,j)=Ef(1,j-1)+x(j)*x(j);
end
gama(1,:)=ones(1,N);
%LSL算法
for m=1:(M-1)
    for n=2:N
        delta(m+1,n)=delta(m+1,n-1)+eb(m,n-1)*ef(m,n)/gama(m,n-1);
        ef(m+1,n)=ef(m,n)-delta(m+1,n)*eb(m,n-1)/Eb(m,n-1);
        eb(m+1,n)=eb(m,n-1)-delta(m+1,n)*ef(m,n)/Ef(m,n);
        Ef(m+1,n)=Ef(m,n)-delta(m+1,n)^2/Eb(m,n-1);
        Eb(m+1,n)=Eb(m,n-1)-delta(m+1,n)^2/Ef(m,n);
        gama(m+1,n-1)=gama(m,n-1)-eb(m,n-1)^2/Eb(m,n-1);
        kf(m,n)=delta(m+1,n)/Ef(m,n);
        kb(m,n)=delta(m+1,n)/Eb(m,n-1);
    end
end
%求全矢量系數
for z=1:N
    a(1,z)=kb(1,z)-kf(1,z)*kb(2,z);
    a(2,z)=kb(2,z);
end
%求出y(n)和d(n)
for n=M:(N-1)
    xl=x(n-1:-1:n-M+1);
    y(n)=a(:,n)'*xl';
    error(n)=x(n)-y(n);   
end
%畫出LSL算法權矢量過渡過程的圖形
a1=a(1,:);
a2=a(2,:);
plot(a1,'r');
hold on;
plot(a2,'b');

x(1)=A(1);
x(2)=A(2)+1.558*A(1);
for i=3:N
    x(i)=A(i)+1.558*x(i-1)-0.81*x(i-2);    
end

d=x;                                    %期望響應d(n)=x(n)
w=zeros(2,N);                           %權矢量w(n)初始化
mu=0.005;                               %步長因子mu
for n=3:N-1                             %LMS算法
    xl=x(n-1:-1:n-2);
    y(n)=w(:,n)'*xl';
    e(n)=d(n)-y(n);
    w(:,n+1)=w(:,n)+2*mu*e(n)*xl';     
end
%畫出LMS算法權矢量過渡過程的圖形
a1=w(1,:);
a2=w(2,:);
plot(a1,'y');
plot(a2,'k');
legend('LSL a1(n)','LSL a2(n)','LMS a1(n)','LMS a2(n)');
plot(z1,'r:');
plot(z2,'r:');
title('LSL和LMS算法比較');
hold off;