ar.m

随机信号处理AR自回归模型

%--------------------------------－－－－－－－－--%
%      随机信号处理上机作业 ： AR模型参数估计       %
% name:  张磊          number:0610120653          %
% 06.12.2                                         %
%-------------------------------－－－－－－－－---%
%                 LMS simulation                  %
%-------------------------------------------------%
a1=0.1;
a2=-0.8;
n=2048;
v=randn(1,n*2);
x=zeros(1,n*2);
for k=1:n*2
    if k==1
        x(k)=randn(1);
    elseif k==2
        x(k)=-a1*x(k-1)+v(k);
    else
        x(k)=-a1*x(k-1)-a2*x(k-2)+v(k);
    end
end

Rxx=zeros(512,512);
X=zeros(512,n);
for i=1:512
    i
    X(i,1:n)=x(1,i:n+i-1);
end
for i=1:50
    for k=1:50
        Rxx(i,k)=sum(X(i,:).*conj(X(k,:)))/n;   %在计算自相关矩阵时，选取512点信号求相关，以时间平均代替
    end                                         %概率平均取得较准确的相关系数估计
end 
%根据LMS准则建立并解YULE-WALKER方程
rxx=zeros(3,3);
rxx(:,:)=Rxx(1:3,1:3);
b=[1 zeros(1,2)]';
Rinv=inv(rxx);
a=Rinv*b;    
% ans:     1.0325    0.1037  -0.8368

%------------------------------------------------------%
%                   RLS估计参数                        %
%------------------------------------------------------%
n=100;          % n 循环次数
W=zeros(3,1); 
W(1)=0;
W(2)=0;
W(3)=1;
R=zeros(3,3);
for k=1:3
    R(k,k)=1;
end
XX=zeros(3,1);
for i=1:n
    if i==1
        XX(1)=0;
        XX(2)=0;
        XX(3)=v(1);
        
        invR=R-R*XX*XX'*R/(1+XX'*R*XX);
        e(i)=x(i)-XX'*W;
        W=W+invR*XX*e(i);
        W(3)=1;
    elseif i==2
        XX(1)=x(1);
        XX(2)=0;
        XX(3)=v(2);
        
        invR=invR-invR*XX*XX'*invR/(1+XX'*invR*XX);
        e(i)=x(i)-XX'*W;
        W=W+invR*XX*e(i);
        W(3)=1;
    else
        XX(1)=x(i-1);
        XX(2)=x(i-2);
        XX(3)=v(i);
        
        invR=invR-invR*XX*XX'*invR/(1+XX'*invR*XX);
        e(i)=x(i)-XX'*W;
        W=W+invR*XX*e(i);
        W(3)=1;     
    end
end
% ans: 加权系数:   
%   -0.1012
%   0.7967
%   1.0000


%-----------------------------------－－－－－－-----%
%             Levison-Durin 算法                    %
%------------------------------------－－－－－－----%

a=zeros(20,20);
sita=zeros(20,1);
FPE=zeros(20,1);
for i=1:20
    if i==1
        a(1,1)=-Rxx(1,2)/Rxx(1,1);
        sita(1)=(1-a(1,1)*(conj(a(1,1))))*Rxx(1,1);
        FPE(i)=sita(i)*(n+i-1)/(n-i+1);
    else
        temp=0;
        for k=0:i-2
            temp=temp+a(i-1,k+1)*Rxx(1,i-k);
        end
        a(i,i)=-(Rxx(1,i+1)+temp)/sita(i-1);
        for k=1:i-1
            a(i,k)=a(i-1,k)+a(i,i)*a(i-1,i-k);
        end
        sita(i)=(1-a(i,i)*conj(a(i,i)))*sita(i-1);
        FPE(i)=sita(i)*(n+i-1)/(n-i+1);  %使用最终预测误差准则FPE 确定AR模型阶数
    end
end
[minFPE index]=min(FPE);
order=index;  % ans: AR模型阶数 = 2  系数： 0.0785   -0.8090