p=input('請輸入AR模型的階數p=');
a=zeros(1,p);
for i = 1:p,                        %根據所輸入的階數輸入參數
    teStr=sprintf('請輸入AR模型的參數a%d=',i);
    a(i)=input(teStr);
end
tempA=zeros(1,p+1);                 %生成最高階數為p的多項式，系數從高次
                                    %到低次依次為ap到a1，常數項為1
tempA(p+1)=1;
for i=1:p,
    tempA(p+1-i)=a(i);
end  
rootp=zeros(p,1);
rootp=roots(tempA)                   %求多項式的所有根
checkA=zeros(p,1);
for i=1:p,
checkA(i,1)=abs(rootp(i,1));       %求多項式的所有根的模
end
checkA
for i=1:p,             %檢驗多項式的所有根的模是否大于1，若全部大于1，繼續執行下面的程序，否則退出
 if check(i,1)<=1; 
    display('s(n)不收斂，錯誤的輸入,請重新輸入');
    return;
    end
end
%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);
s=zeros(1,L);
for i = 1:L,                                  %計算出s(n)
    te=0;
    for m = 1:p,
        if i-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)=deltaW;
A=Rss^(-1)*epsilon;
for i = 1:p,     
    sprintf('a%d=%f\n',i,A(i+1))
end
deltaW_e=0;
for i = 1:p,
    deltaW_e=deltaW_e+a(i)*Rss(1,i+1);
end
deltaW_e=deltaW_e+Rss(1,1);
sprintf('deltaW_e=%f',deltaW_e)