单点平稳模拟中，如何通过matlab绘制出多个样本函数的集合自相关函数图像？

问题遇到的现象和发生背景

在平稳单点模拟中，如何通过matlab画出多个样本函数的集合自相关函数图像？下面的模拟公式怎么用更快的方式来编程实现，由于本人才刚开始学习matlab，编的代码感觉效率太低。

自己编的matlab代码

clc
clear all
o=1;b=1;
wu=4*pi;
Nw=500;
Ns=128;
N=256;
dw=wu/Ns;
T0=50;
dt=0.25;
w=0:dw:4*pi-dw;
t=0:dt:T0-dt;
T=0:dt:T0-dt;
RT=o.^2.*b.^4.*(b.^2-3*T.^2)./((b.^2+T.^2).^3);
phi=rand(Nw,N)*2*pi;
wn=zeros(Nw,N)
for n=1:1:N;
wn(:,n)=n.*dw;
end
Sff=0.25.*o.^2.*b.^3.*wn.^2.*exp(-b.*abs(wn));
An=(2.*Sff.*dw).^0.5;
f=zeros(Nw,N);
M=zeros(1,N);
for j=1:1:Nw;
for n=1:1:N
f1=sqrt(2).*An(j,n).*cos(wn(j,n).*t+phi(j,n));
M=f1+M
end
f(j,:)=M;
end
plot(t(1:N),f(1,1:N))
maxlag=N;
estacor_1=xcorr(f(1,:),maxlag,'biased')
tao=-maxlag/4:1/4:maxlag/4;
hold on
plot(tao,(estacor_1)-(max(estacor_1)-max(RT)),'r--','LineWidth',2.8);
hold on
plot(T,RT,'k','linewidth',2.5);
legend('f','RT*','RT')
axis([0 64 -3 3])

循环赋值那一行代码效率很低，很费时；画出的图像与论文要求的RT图像没有完全吻合

后面的画图自相关函数是我找师兄问的，但是没怎么看懂希望能帮我解释一下，或者提供一个画出集合自相关函数图像的方法。

我想要达到的结果

1.集合自相关函数图像的绘制方法；
2.优化代码；

clear all
clc
T=0.05;             %采样周期
fs=1/T;             %采样频率
t=-200:T:200;
fi=pi*(0:0.1:2);    %相位
h=2*t+5;
figure,plot(t,h),title('线性系统')
for i=1:length(fi)
    x(i,:)=cos(t+fi(i));
    y(i,:)=conv(x(i,:),h);
end
%--------------------------------正余弦信号经过线性系统后的包络
tt=-200:0.05:600;                     %卷积后时域横轴坐标
x1=sin(t);
y1=conv(x1,h);
x2=cos(t);
y2=conv(x2,h);
figure,plot(tt,y1),xlim([-240,650])
title('正弦信号经过线性系统'),xlabel('t'),ylabel('y(t)')
figure,plot(tt,y2),xlim([-240,650])
title('余弦信号经过线性系统'),xlabel('t'),ylabel('y(t)')
% ------------------正弦信号相位不同，经过线性系统后输出的包络不同
for i=1:length(fi)
    figure
    plot(tt,y(i,:)),xlim([-240,650])
end

解题思路：
通过调用一次 plot，多个 x-y 对组参数会创建多幅图形。MATLAB对每条线使用不同的颜色。【利用枚举法】
参考链接【平稳随机过程的样本函数（随机信号）通过线性系统的仿真】，期望对你有所帮助：https://www.freesion.com/article/3158432798/
【仅供参考】