%%信號源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