p=input('請輸入AR模型的階數p=');
a=zeros(1,p);
for i = 1:p,                        %根據所輸入的階數輸入參數
    teStr=sprintf('請輸入AR模型的參數a%d=',i);
    a(i)=input(teStr);
end
deltaW=input('請輸入白噪聲方差:');
L=input('請輸入信號s(n)樣本個數L:');
w=sqrt(deltaW)*randn(1,L);          %產生方差為deltaW的高斯白噪聲
s=zeros(1,L);
s(1)=w(1);                          %根據式1－20生成s(n)
for i = 1:L,          
    te=0;
    for m = 1:p,
        if i-m>0,                   %凡是n-m<＝0的一律s(n-m)=0      
           te=te+a(m)*s(i-m);
        end
    end
    s(i)=w(i)-te;
end
Rss = zeros(p+1,p+1);               %Rss為s(n)的自相關
for i = 1:p+1,                      %i為行數
    for j = 1:p+1,                  %j為列數
        teR=0;   
        te=abs(i-j);                %括號中的值是行數減去列數
        for k = 1:L-te,             %由式1-18得到Rss 
            teR=teR+s(k)*s(k+te);
        end
        teR=teR/(L-te);
        Rss(i,j)=teR;               %賦值給Rss
    end
end
A=zeros(p+1,1);
epsilon=zeros(p+1,1);
epsilon(1)=1;
for i=2:(p+1),
    epsilon(i)=0;
end
A=Rss^(-1)*deltaW*epsilon;                 %由式1-21可得到
for i = 1:p,     
    sprintf('估計的a%d=%f\n',i,A(i+1))
end
deltaW_e=Rss(1,1);
for i = 1:p,
    deltaW_e=deltaW_e+a(i)*Rss(1,i+1);
end
sprintf('估計的deltaW_e=%f',deltaW_e)