%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%題目：????????????????????????????????????????????????????????????????????%
%已知AR（4）過程：???????????????????????????????????????????????????????????%
%x（n）=2.760x（n-1）-3.809x（n-2）+2.654x（n-3）-0.924x（n-4）+w（n）????????????????%
%?試用Burg算法利用Matlab估計模型參數、階數和功率譜。?????????????????????????%
%實驗報告要求：實驗仿真100次，但不要求給出所有100次結果，?????????????????????%
%只給出估計結果的方差和偏差（用估計結果均值與真值的差來估計偏差）。????????????%
%實驗報告要求寫清楚基本原理，提交程序并要有關鍵注釋。?????????????????????????%
%電子通信工程?研1205?于飛?2012020148????????????????????????????????????????%
%2012年10月11日%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear?allclose?allclc;

sum=zeros（5001）;%psd估計值寄存器
for?o=1:1:100
????%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
????%產生均值為零、方差為1、服從高斯（正態）分布的白噪聲信號u（n）%
????%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
????w=randn（150000）;
????w=w-mean（w）;
????power_u=var（w）;
????m_u=mean（w）;
????%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
????%使u（n）經過系統?x（n）=2.760x（n-1）-3.809x（n-2）+2.654x（n-3）-0.924x（n-4）+w（n）%
????%得到所要研究的隨機過程x（n）%
????%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
????c=1;

????%d=[1-1.3521.338-0.6620.240];%寬帶過程
????d=[1-2.7603.809-2.6540.924];%窄帶過程
????y=filter（cdw）;
????%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
????%burg算法計算參數、階數、功率譜?????????????????????????????????????????????%
????%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
????N=256;
????x=y（25000:25000+N-1）;%截取長度為256的樣本
????f=x;%初始條件，前向預測值
????g=x;%初始條件，后向預測值
????p=（1/N）*x*x‘;%初始誤差功率
????a=zeros（110）;%產生用來估計系數的向量a
????b=zeros（110）;%產生用來估計系數的向量b
????a（1）=1;%初始條件，初始系數a0=1;

????for?m=1:1:10%循環遞推，最多10階
????????num=-2*f（m+1:N）*g（m:N-1）‘;%計算反饋系數k的分子
????????den=f（m+1:N）*f（m+1:N）‘+g（m:N-1）*g（m:N-1）‘;%計算反射系數k的分母
????????k=num/den;%計算反射系數
????%%%%%%%%%%%%%%%%%%%%%%-----------計算系數a------------%%%%%%%%%%%%%%%%%%%%%%
????????for?n=1:1:m-1
???????????b（n）=a（m-n）;
????????end;
????????if?m==1
????????????a（1）=k;
????????else
????????????b（m）=1;
????????end;

????????a=a+k*b;
????%%%%%%%%%%%%%%%%%%%%%%------------------------------%%%%%%%%%%%%%%%%%%%
????%%%%%%%%%%%%%%%%%-----------計算前、后項預測值-----------%%%%%%%%%%%%%%%
????????ft=f;
????????gt=g;
????????f=ft+k*[0gt（1:N-1）];
????????g=k*ft+[0gt（1:N-1）];
????%%%%%%%%%%%%%%%%%%%%%%------------------------------%%%%%%%%%%%%%%%%%%%
????????pn=（1-k^2）*p;%計算誤差功率
????????if?p-pn?????????????break;
????????end;
????????p=pn;
????end;
????figure（1）
????a=[1a（1:m）];
????[HW]=freqz（1a500）;
????psd=power_u*（abs（H）.^2）;
????sum?=sum?+psd;
????%plot（（1/（2*pi））*W10*log10（psd））;%用估計的參數估計功率譜并繪出藍色
????hold?on;
????if?o==1
????????[HW]=freqz（cd500）;
????????psdn=power_u*（abs（H）.^2）;
????????plot（（

?屬性????????????大小?????日期????時間???名稱
-----------?---------??----------?-----??----
?????目錄???????????0??2012-11-10?10:51??burgyf\
?????文件????????3493??2012-10-14?14:24??burgyf\ex2.asv
?????文件????????3537??2012-10-14?14:30??burgyf\ex2.m
?????文件???????69120??2012-11-08?22:25??burgyf\實驗報告.doc