AR 谱估计 matlab版

%%信号源1
clear all;
close all;
clc;

N=512;%采样点数
MU=0;%噪声均值
D=0.101;%噪声方差
times=100
for i=1:times
v=normrnd（MUsqrt（D）N1）;%随机产生均值为MU方差为D的高斯白噪声
fs=1;%采样频率为1hz
t=1/fs;%采样周期
f1=0.05;
f2=0.40;
f3=0.42;

z=0;
a1=-0.85;
for i=1:N
    z=-1*a1*z+v（i）;
    x（i）=cos（2*pi*f1*i）+cos（2*pi*f2*i）+cos（2*pi*f3*i）+0.5*z;
end


 p=12;


Rx=xcorr（x‘biase‘）;
Rxx=toeplitz（Rx（N:N+p-1））;
r=-Rx（N+1:N+p）;
a=inv（Rxx）*r‘;

w=linspace（0pi256）;
f=w/（2*pi）;
e=zeros（1length（w））;
for i=1:p;
    e=a（i）*exp（-j*w*i）+e;
end
H=abs（e+1）.^2;
P=D./H;
P=10*log10（P）;
subplot（222）
plot（fP）
xlabel（‘频率HZ‘）;
ylabel（‘功率db‘）;
title（‘信号源1的估计功率谱图‘）;
 hold on;

end
%%

%%信号源二
clear all;
close all;
clc;

N=64;
MU=0;
D=1;
times=100
for i=1:times
v=normrnd（MUsqrt（D）N1）;%随机产生均值为MU方差为D的高斯白噪声
fs=1;%采样频率为1hz
t=1/fs;%采样周期
% f1=0.05;
% f2=0.40;
% f3=0.42;
z=0;

a1=-1.35;
a2=1.34;
a3=-0.66;
a4=0.24;
x=zeros（N+41）;

for i=5:N+4
x（i）=-a1*x（i-1）-a2*x（i-2）-a3*x（i-3）-a4*x（i-4）+v（i-4）;
 
end
x（[1:4]）=[];
%  p=computeAR_P（xN）;
 p=7;
Rx=xcorr（x‘biase‘）;
Rxx=toeplitz（Rx