%设定相关的常量
sigma = 1;
a = 0.8;
%设置产生的序列元素个数为500个
N = 500;
%将C视作对角阵D，下三角矩阵L和上三角矩阵L'的相加
D = eye(N);
L = zeros(N);
for i = 2:500
    for j = 1:i-1
        L(i,j) = a.^(i-j);
    end
end
C = sigma^2/(1-a^2).*(D + L' +L);
%将协方差矩阵C进行分解
A = (chol(C))';

%获得标准正态随机矢量
U = normrnd(0,1,500,1);
%获得相关正态随机序列并打印出来
X = A*U;
plot(X);
title('相关正态随机序列')
print -dpng 报告/figures/Correlated_normal_random_sequence.png

%实际相关函数
R = xcorr(X,'coeff');
plot(R,'b','LineWidth',1.0);
hold on;
%理论相关函数
m = -500 : 500;
RX = 1/(1-0.64).*(0.8).^(abs(m));
plot(RX,'r','LineWidth',1.0);
title('相关函数')
legend('实际曲线','理论曲线')
print -dpng 报告/figures/correlation_function.png
hold off

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