clear
clc
%系统的初始输入，为平稳白噪声的一个值
W = 0+2.*randn(1,500);
X = [W(1)];
%模拟AR过程产生500个观测点
for n = 2:500
    X = [X, -0.8*X(n-1)+W(n)];
end
%画出该样本函数的波形
n = 1:500;
plot(n,X);
title('500个观测点的波形')
print -dpng 报告/figures/500_observation_point_waveform.png
%计算序列的均值和方差
x_mean = mean(X)
x_variance = var(X)

%理论自相关函数
a=-0.8;
m = -50:0.1:50;
R = 4.*a.^abs(m)/(1-a.^2);
plot(m,R);
title('理论自相关函数')
print -dpng 报告/figures/Theoretical_autocorrelation_function.png

%理论功率谱
a = -0.8;
omega = -6*pi:0.1:6*pi;
G = 4./(1+a.^2-2*a.*cos(omega));
plot(omega,G);
title('理论功率谱')
print -dpng 报告/figures/Theoretical_power_spectrum.png

%实际自相关函数
r = xcorr(X,'biased');
plot(r);
title('实际自相关函数')
print -dpng 报告/figures/autocorrelation_function_for_three.png

%实际功率谱
 window = hann(length(X));
[Pxx,w] = periodogram(X,window);
plot(w,Pxx);
title('实际功率谱密度')
print -dpng 报告/figures/power_spectrum_for_three.png