S=sqrt(0.05)*randn(1000,1);     %產生一個白噪聲序列S(n)
H=[3; -4];                      %初始濾波器系數H(0)
T=[0 0];
T2=[0 0];
e=0;
e2=0;
for n=0:1:300                   %LMS迭代
    y=sin(pi/8*(n+1)+pi/10)+S(n+1,:);
                                %計算每步參考信號y(n+1)
    e=y-H'*[sqrt(2)*sin(pi/8*(n+1)); sqrt(2)*sin(pi/8*n)];
                                %計算每步誤差e(n+1)
    e2=e^2;                     %計算每步J(n+1)=e^2(n+1)
    H=H+0.4*e*[sqrt(2)*sin(pi/8*(n+1)); sqrt(2)*sin(pi/8*n)];
                                %計算每步濾波器系數H(n+1)
    T=[T; [n+1,e]];
    T2=[T2; [n+1,e2]];
end

%繪制e(n)波形
figure(1);
hold on;
axis([0,300,-2,2]);
T=T';
plot(T(1,:),T(2,:));

%繪制J(n)=e^2(n)隨時間n的變化曲線
figure(2);
hold on;
axis([0,300,0,4]);
T2=T2';
plot(T2(1,:),T2(2,:));

%繪制J(n)=e^2(n)取100次實驗結果平均值的變化曲線
figure(3);
hold on;
axis([0,300,0,4]);
A=[1:301; linspace(0,0,301)];
A=A';
for i=0:99                      %100次LMS迭代實驗
    S=sqrt(0.05)*randn(1000,1);
    H=[3; -4];
    T=[0 0];
    T2=[0 0];
    e=0;
    e2=0;
    for n=0:1:300
        y=sin(pi/8*(n+1)+pi/10)+S(n+1,:);
        e=y-H'*[sqrt(2)*sin(pi/8*(n+1)); sqrt(2)*sin(pi/8*n)];
        e2=e^2;
        H=H+0.4*e*[sqrt(2)*sin(pi/8*(n+1)); sqrt(2)*sin(pi/8*n)];
        A(n+1,:)=A(n+1,:)+[0,e2]/100;
                                %所有J(n)實驗結果取平均值
    end
end
A=A';
plot(A(1,:),A(2,:));            %繪制J(n)平均值的變化曲線