g=100;              %統計仿真次數為g
N=256;              %輸入信號抽樣點數
k=128;               %時域抽頭LMS算法濾波器階數
pp=zeros(g,N-k);    %將每次獨立循環的誤差結果存于矩陣pp中，以便后面對其平均
u=0.01;             %濾波器收斂因子
for q=1:g
    t=1:N;
    a=1;
    s=a*sin(0.05*pi*t);    %輸入單信號s
    figure(1);
    subplot(311)
    plot(t,real(s));       %信號s時域波形
    title('信號s時域波形');
    xlabel('n');
    ylabel('s');
    axis([0,N,-a-1,a+1]);
    xn=awgn(s,5);           %加入均值為零的高斯白噪聲，信噪比為3dB
    %設置初值
    y=zeros(1,N);           %輸出信號y
    y(1:k)=xn(1:k);         %將輸入信號xn的前k個值作為輸出y的前k個值
    w=zeros(1,k);           %設置抽頭加權初值
    e=zeros(1,N);           %誤差信號
   %用LMS算法迭代濾波
    for i=(k+1):N
        XN=xn((i-k+1):(i));
        y(i)=w*XN';
        e(i)=s(i)-y(i);
        w=w+u*e(i)*XN;
    end
    pp(q,:)=(e((k+1):N)).^2;
end
subplot(312)
plot(t,real(xn));                   %信號s時域波形
title('信號s加噪聲后的時域波形');
subplot(313)
plot(t,real(y));                    %信號s時域波形          
title('自適應濾波后的輸出時域波形');

for b=1:N-k
    bi(b)=sum(pp(:,b))/g;           %求誤差統計平均
end
figure(2);                          %算法收斂曲線
t=1:N-k;
plot(t,bi,'r');
hold on                             %將每次循環的圖形顯示結果保存下來