%基于RLS算法的自適應線性預測 
clc; 
clear all; 
N=300; 
M=100;%計算的次數 
w1=zeros(N,M);w2=zeros(N,M);I=eye(2);e1=zeros(N,M); 
 
for k=1:M 
    %產生白噪聲 
    Pv=0.008;%定義白噪聲方差 
    a1=-0.195;a2=0.95;o=0.02;r=0.95; 
    m=5000;%產生5000個隨機數 
    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(2)-a1*v(1);%x(1)=v(1)-a1*v(0) 
    w=zeros(2,N); 
    w(:,1)=[0 0]';%w(0)=[0 0]'; 
    X=zeros(2,N); 
    X(:,2)=[v(1) 0]';%X(0)=[0 0]';X(1)=[v(0) 0]' 
    C=zeros(2,2*N); 
    C(:,1:2)=1/o.*I;%C(0)=1/o*I 
    e=zeros(1,N)';%定義誤差向量 
    u=zeros(1,N); 
    g=zeros(2,N); 
 
    %根據RLS算法進行遞推 
 
    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)]'; 
        u(n)=X(:,n+1)'*C(:,2*n-1:2*n)*X(:,n+1); 
        g(:,n)=(C(:,2*n-1:2*n)*X(:,n+1))./(r+u(n)); 
        w(:,n+1)=w(:,n)+g(:,n)*(x(n+1)-X(:,n+1)'*w(:,n)); 
        C(:,2*n+1:2*(n+1))=1/r.*(C(:,2*n-1:2*n)-g(:,n)*X(:,n+1)'*C(:,2*n-1:2*n)); 
        e(n)=x(n+1)-X(:,n+1)'*w(:,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)的單次變化曲線(線性預測，RLS)') 
subplot(222) 
plot(n,wa1(n,1)./M,n,wa2(n,1)./M);%作出100次計算權矢量的平均變化曲線 
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的平均變化曲線');