%AR（3）產(chǎn)生數(shù)據(jù)，用3階的LD算法計算，輸出估計的功率譜
clear;
p=100;??????????%AR計算的階數(shù)

u=wgn（100011）;???%u（n）是一個均值為零，功率為一的白噪聲
x（1）=u（1）;
x（2）=u（2）;
x（3）=u（3）;
for?kk=4:1000
???x（kk）=14/24*x（kk-1）+9/24*x（kk-2）-1/24*x（kk-3）+u（kk）;
end???????????%產(chǎn)生指定的AR（3）模型的x（n）
figure?
pburg（x3）
R?=?zeros（1?p+1）;
???sample_num?=?length（x）;
???for?n?=?1:p+2
???????R（n）?=?mean（?x（1:sample_num-n+1）?.*?x（n:sample_num）?）;%自相關(guān)函數(shù)
???end
%for??n=?2:p+2
???%?R（n）=R（n）/R（1）;
%end
????ssigma?=?zeros（1?p）;???%?噪聲方差矩陣
????D?=?zeros（1?p）;????????%?Dk矩陣
????gama?=?zeros（1?p）;?????%?gama矩陣
????a?=?zeros（p?p）;???????
???
????a0?=?1;
????D0?=?a0?*?R（1+1）;
????ssigma0?=?R（1）;
????gama（1）?=?D0?/?（ssigma0?+?eps）;
????a（11）?=?-gama（1）;
????
????for?k?=?1:p
????????%公式中的?R（k）?is?R（k+1）?in?this?program
????????ssigma（k）?=?R（1）?+?sum（?a（k?1:k）?.*?R（2:k+1）?）;
????????D（k）?=?a0*R（k+2）?+?sum（?a（k?1:k）?.*?fliplr（R（2?:?k+1