AR模型功率谱估计matlab仿真

含注释，参数模型功率谱估计AR模型自相关法仿真，原理来自《数字信号处理理论、算法与实现》第三版：P545-P547

%%%%%%%%%%%%%%%%%%%现代谱估计///自相关法%%%%%%%
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）;         %相关系数
rx=R_X（1N+1:2*N-1）;
g=zeros（1p）;         %误差功率
K=zeros（1p）;         %反射系数
a=zeros（pp）;         %模型系数（方阵）
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）; %系数
        end
end
%迭代结束，总结结果
Ar=a（p1:p）;                  %AR模型系数（第p行所有的系数）
Ar=[1Ar];                    %1*P+1矩阵
G=g（1p）;                     %最小误差功率
 F_Ar=abs（fft（ArN））;         %系数的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模型自相关法阶数数据长为1024，P=‘ int2str（p）] ）; 

 [y3fy]=pyulear（ypNfs）;
 subplot（212）
 plot（fy10*log（y3））
 xlabel（‘频率‘）;ylabel（‘功率谱P（dB）‘）;
title（‘matlab自带函数pyulear‘）; 

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     文件        1490  2019-11-07 21:25  zixiangguan.m