num=[1 7 24 24]; den=[1 10 35 50 24]; G=tf(num,den);
[B,A]=spec_fac(3*num,den); G1=tf(B,A);
[y0,t]=impulse(G1,6); % 求取系統的脈沖響應數據
t=[-t(end:-1:2); t]; y0=[y0(end:-1:2); y0]/max(y0); % 鏡像延拓并歸一化

G=tf([1,7,24,24],[1,10,35,50,24]); T=0.02;
[F,G0,D,C]=sc2d(G,3,T);
n_point=30000; r=randn(n_point+4,1); r=r-mean(r);
y=zeros(n_point,1); x=zeros(4,1); d0=0;
for i=1:n_point
   x=F*x+G0*d0+D*r(i:i+3); y(i)=C*x;
end % 以上語句可以得出仿真結果
[Cyy,f]=crosscorr(y,y,300); % 由互相關函數求自相關函數，自動延拓
f=f*T; plot(f,Cyy,t,y0,':') % 繪制曲線，理論值用虛線表示