%基于LMS算法的自適應線性預測

clc;clear all;
N=1000;
M=100;%計算的次數
w1=zeros(N,M);w2=zeros(N,M);e1=zeros(N,M);

for k=1:M
    %產生白噪聲
    Pv=0.0965;%定義白噪聲方差
    a1=-0.195;a2=0.95;u=0.01;
    m=50000;%產生50000個隨機數
    v=randn(1,m);
    v=v*sqrt(Pv);%產生均值為0，方差為Pv的白噪聲
    
    m=1:N;
    v=v(1:N);%取出前1000個
    %plot(m,v);title('均值為0，方差為0.0965的白噪聲');ylabel('v(n)');xlabel('n');
    v=v';

    %向量初使化
    x=zeros(1,N)';
    x(1)=v(1);%x(0)=v(0)
    x(2)=v(1)-a1*v(1);%x(1)=v(1)-a1*v(0)
    w=zeros(2,N);%w(0)=[0 0]';
    X=zeros(2,N);
    X(:,2)=[v(1) 0]';%X(0)=[v(0),0]';X(1)=[v(0) 0]'
    e=zeros(1,N)';%定義誤差向量

    %根據LMS算法進行遞推

    for n=1:N-2
        x(n+2)=v(n+2)-a1*x(n+1)-a2*x(n);
        X(:,n+2)=[x(n+1) x(n)]';
        e(n)=x(n+1)-X(:,n+1)'*w(:,n);
        w(:,n+1)=w(:,n)+u*X(:,n+1)*e(n);       
        
        w1(:,k)=w(1,:)'; w2(:,k)=w(2,:)';%將每次計算得到的權矢量值儲存 
        e1(:,k)=e(:,1);%將每次計算得到的誤差儲存
     end
     
 end

%求權矢量和誤差的M次的平均值         
wa1=zeros(N,1);wa2=zeros(N,1);en=zeros(N,1); 
for k=1:M
    wa1(:,1)=wa1(:,1)+w1(:,k);
    wa2(:,1)=wa2(:,1)+w2(:,k);
    en(:,1)=en(:,1)+e1(:,k);
end

n=1:N;
subplot(221)
plot(n,w(1,n),n,w(2,n));%作出單次計算權矢量的變化曲線
xlabel('n');ylabel('w(n)');title('w1(n)和w2(n)的單次變化曲線(線性預測，LMS)');
subplot(222)
plot(n,wa1(n,1)./M,n,wa2(n,1)./M);%作出M次計算權矢量的平均變化曲線
xlabel('n');ylabel('w(n)');title('w1(n)和w2(n)100次平均變化曲線');
subplot(223)
plot(n,(e(n,1)).^2);%作出單次計算e^2的變化曲線
xlabel('n');ylabel('e^2');title('單次計算e^2的變化曲線');
subplot(224)
plot(n,(en(n,1)./M).^2);%作出M次計算e^2的平均變化曲線
xlabel('n');ylabel('e^2');title('100次計算e^2的平均變化曲線');