function?ARMODEL（）
%現(xiàn)代數(shù)字信號處理作業(yè)
%采用MATLAB仿真實(shí)現(xiàn)AR模型，進(jìn)行譜估計(jì)
%AR模型的理論公式：?x（n）?+?a1*x（n-1）?+?a2*x（n-2）?+?……?+ap*x（n-p）?=?w（n）
%待估計(jì)數(shù)據(jù)樣本：
%目前輸入隨機(jī)過程為疊加了正態(tài)分布噪聲的正弦波

Fs?=?1000;??
t?=?0:1/Fs:15;
N?=?size（t2）??????????????????????%數(shù)據(jù)樣值點(diǎn)數(shù)
randn（‘state‘0）;
x?=?cos（2*pi*t*200）?+?randn（1N）;??%?200Hz?cosine?plus?noise

%計(jì)算N個(gè)取樣數(shù)據(jù)的取樣數(shù)據(jù)自相關(guān)函數(shù)
rxx?=?zeros（1N）;???%保存取樣數(shù)據(jù)自相關(guān)函數(shù)的變量
for?m?=?0:N-1
????sum?=?0;
????for?n=?1:N-m
????????temp1?=?x（n）*x（m+n）;
????????sum?=?sum?+?temp1;
????end
????rxx（m+1）?=?sum/N;
end

%采用Levison-Durbin算法求解AR模型的Yule-Walker模型
%需要確定AR模型理論公式中的參數(shù)：白噪聲w（n）的方差、方程系數(shù)a1……ap（這里包括了模型的階次）
PMAX?=?100;????????????????%設(shè)定AR模型最高階次
atemp1?=?zeros（1PMAX+1）;?%保存方程系數(shù)的中間變量
atemp2?=?zeros（1PMAX+1）;?%保存方程系數(shù)的中間變量
deviationtemp1?=?zeros;?%保存白噪聲w（n）方差的中間變量
deviationtemp2?=?zeros;?%保存白噪聲w（n）方差的中間變量

%AR（1）模型：x（n）?+?a1*x（n-1）?=?w（n）
%其Yule-Walker方程：?R（0）*1?+?R（1）*a1?=?deviation1;
%????????????????????R（1）*1?+?R（0）*a1?=?0;
%求解方程確定a1、deviation1
atemp1（1）?=?1;
atemp1（2）?=?-rxx（2）/rxx（1）;
atemp2?=?atemp1;
deviationtemp1?=?（????rxx（1）*rxx（1）???-???rxx（2）*rxx（2）???）/rxx（1）;
deviationtemp2?=?deviationtemp1;

%利用Levison-Durbin迭代算法計(jì)算AR模型參數(shù)
%根據(jù)FPE準(zhǔn)則、AIC準(zhǔn)則和BIC準(zhǔn)則確定AR模型的階次
%atemp1、deviation1保存第k次的運(yùn)算結(jié)果
%atemp2、deviation2保存第k+1次的運(yùn)算結(jié)果
FPE（1）?=?deviationtemp1*（N+2）/N;
AIC（1）?=?log（deviationtemp1）?+?2/N;
BIC（1）?=?log（deviationtemp1）?+?log（N）/N;
veriance（1）?=?deviationtemp1;

criteria?=?3?

for?P?=?2:PMAX
????sum1?=?0;
????sum2?=?0;
????for?i?=?2:（P+1）
????????sum1?=?atemp1（i）*rxx（i）?+?sum1;
????end
????
????for?i?=?1:（P+1）
????????sum2?=?atemp1（i）*rxx（P+2-i）?+?sum2;
????end
????
????deviationtemp1?=?rxx（1）?+?sum1;
????dk?=?sum2;
????ref（P）?=?dk/deviationtemp1;
????deviationtemp2?=?（?1?-?ref（P）*ref（P）?）*deviationtemp1;
????
????for?i?=?2:（P+1）
????atemp2（i）?=?atemp1（i）?-?ref（P）*atemp1（P+2-i）;
????end
????%計(jì)算AR（P）模型參數(shù)
????atemp1?=?atemp2;
????veriance（P）?=?deviationtemp2;
????
????%利用階次判斷準(zhǔn)則
????FPE（P）?=?veriance（P）*（1+P/N）/（1-P/N）;
????AIC（P）?=?log（veriance（P））?+?2*P/N;
????BIC（P）?=?log（veriance（P））?+?P*log（N）/N;
????
	if?（criteria?==?1）?&?（FPE（P-1）????????%利用最終預(yù)測誤差（FPE）準(zhǔn)則作為AR模型階次的判斷準(zhǔn)則??????
????????a?=?atemp2（1:（P））;
????????P?=?P?-?1;
????????deviation?=?veriance（P-1）;
????????break;??
	end
	if?（criteria?==?2）?&?（AIC（P-1）?????????%利用AIC準(zhǔn)則作為AR模型階次的判斷準(zhǔn)則
????????a?=?atemp2（1:（P））;
????????P?=?P?-?1;
????????deviation?=?veriance（P-1）;
????????break;
	end
	if（criteria?==?3）?&?（BIC（P-1）????????%利用BIC準(zhǔn)則作為AR模型階次的判斷準(zhǔn)則
?????????a?=?atemp2（1:（P））;
?????????P?=?P?-?1;
?????????deviation?=?veriance（P-1）;
?????????break;?
????end
end

%計(jì)算功率譜估計(jì)值
a?=?atemp2（1:（P+1））;
deviation?=?veriance（P）;
freqsample?=?10;
NF?=?Fs/freqsample;
freq?=?[freqsample:freqsample:Fs];
delta_w?=?2*pi/NF;
theta?=?[delta_w:delta_w:2*pi];?
for?k?=?1:NF
????sum?=?0;
????for?i=?2:P+1
????????theta1?=?-1*theta（k）*（i-1）;
????????sum?=?a（i）*（????cos（theta1）+j*sin（theta1）???）?+?sum;
????end

?屬性????????????大小?????日期????時(shí)間???名稱
-----------?---------??----------?-----??----

?????文件???????3735??2004-06-21?18:40??ARMODEL.m

?????文件???????1116??2005-11-13?09:24??zuoye.m

?????文件?????219136??2005-11-23?12:02??隨機(jī)序列作業(yè).doc

-----------?---------??----------?-----??----

???????????????223987????????????????????3