%%%%%%%%%%%%%%%%%%%現(xiàn)代譜估計///自相關(guān)法%%%%%%%
clcclear
N=1024;
n=0:N-1;
fs=1;
x=cos（0.35*pi*n）+cos（0.4*pi*n）+0.25*cos（0.8*pi*n）;
a=0.5;
noise?=?a*randn（size（x））;??%加高斯白噪聲的信號
y?=?x+?noise;
p=input（‘請輸入AR模型的階次（大于2）：‘）;
R_X=xcorr（yN-1‘biased‘）;%xcorr（y‘biased‘）;
rx0=R_X（1N）;?????????%相關(guān)系數(shù)
rx=R_X（1N+1:2*N-1）;
g=zeros（1p）;?????????%誤差功率
K=zeros（1p）;?????????%反射系數(shù)
a=zeros（pp）;?????????%模型系數(shù)（方陣）
g0=rx0;???????????????%初始條件最小誤差功率
K（11）=-rx（11）/rx0;?
a（11）=K（11）;
g（11）=g0*（1-abs（a（11）^2））;
for?j=2:p
?????????K（1j）=-（rx（1j）+sum（a（j-11:j-1）.*rx（1j-1:-1:1）））/g（1j-1）;
?????????a（jj）=K（1j）;
?????????g（1j）=g（1j-1）*（1-K（1j）^2）;?????%誤差功率
????????for?i=1:j-1
????????????a（ji）=a（j-1i）+K（1j）*a（j-1j-i）;?%系數(shù)
????????end
end
%迭代結(jié)束，總結(jié)結(jié)果
Ar=a（p1:p）;??????????????????%AR模型系數(shù)（第p行所有的系數(shù)）
Ar=[1Ar];????????????????????%1*P+1矩陣
G=g（1p）;?????????????????????%最小誤差功率
?F_Ar=abs（fft（ArN））;?????????%系數(shù)的fft
?p_ar=F_Ar.^2;
?P_ar（1:）=G./p_ar（1:）;
?P_ar=10*log10（P_ar）;
?y2=P_ar（1:N/2）;
?f=0:1/N:0.5-1/N;
%??y2=[P_ar（N/2+1:N）P_ar（1:N/2）];
%??f=-0.5:1/N:0.5-1/N;
?%nn=2*pi*f/pi/2;?
?subplot（211）
?plot（fy2）
xlabel（‘頻率‘）;ylabel（‘功率譜P（dB）‘）;
title（[‘AR模型自相關(guān)法階數(shù)數(shù)據(jù)長為1024，P=‘?int2str（p）]?）;?

?[y3fy]=pyulear（ypNfs）;
?subplot（212）
?plot（fy10*log（y3））
?xlabel（‘頻率‘）;ylabel（‘功率譜P（dB）‘）;
title（‘matlab自帶函數(shù)pyulear‘）;?

?屬性????????????大小?????日期????時間???名稱
-----------?---------??----------?-----??----
?????文件????????1490??2019-11-07?21:25??zixiangguan.m