Объединение графиков - Matlab - Ответ 3129568

N=10000                       % число реализаций
del=1/2                       % интервал дискретизации
T=2.5                         % длительность реализаций 
a=2
b=4/2*pi
t=0:del:T;
r=2*exp(-t)-exp(-2*t)         % функция корреляции
n=length(t)                   % число отсчётов 
 for i=1:n
    for j=1:n
        b(i,j)=r(abs(i-j)+1);
    end
end
BB=b                          % исходная корреляционная матрица
D=inv(BB)                     % матрица точности
D(1,6)=0
D(6,1)=0
B=inv(D)                      % аппроксимирующая корреляционная матрица
 
% корреляционные матрицы первых значений
 
B3=[B(1,1)    B(1,2)
        B(2,1)    B(2,2)]
B4=[B(1,1)    B(1,2)    B(1,3)
        B(2,1)    B(2,2)    B(2,3)
        B(3,1)    B(3,2)    B(3,3)]
B5=[B(1,1)    B(1,2)    B(1,3)    B(1,4)
        B(2,1)    B(2,2)    B(2,3)    B(2,4)
        B(3,1)    B(3,2)    B(3,3)    B(3,4)
        B(4,1)    B(4,2)    B(4,3)    B(4,4)]
 
x(1)=randn*sqrt(B(1,1))                         % первое значение последовательности
 
c=B(2,1)/B(1,1)
si=sqrt(B(2,2)-c*B(1,2))
x(2)=c*x(1)+randn*si                            % второе значение последовательности
 
C=[B(3,1) B(3,2)]*inv(B3)
si=sqrt(B(3,3)-C*[B(1,3);B(2,3)])
x(3)=C(1)*x(1)+C(2)*x(2)+randn*si       % третье значение последовательности
 
C=[B(4,1) B(4,2) B(4,3)]*inv(B4)
si=sqrt(B(4,4)-C*[B(1,4);B(2,4);B(3,4)])
x(4)=C(1)*x(1)+C(2)*x(2)+C(3)*x(3)+randn*si           % четвёртое значение
 
%----------------------------------------------------------
                      % установившийся режим генерирования
 
C=[B(5,1) B(5,2) B(5,3) B(5,4)]*inv(B5)               % матрица связи
s=sqrt(B(5,5)-C*[B(1,5);B(2,5);B(3,5);B(4,5)])        % с.к.о. 
X=[x(1);x(2);x(3);x(4)]
for i=5:N
    x(i)=C*X+s*randn;                % следующее значение
    for j=1:4
        X(j)=x(i-4+j);               % сдвиг предыдущих четырёх
    end
end
plot(x)                              % марковская последовательность, 
pause
 
 for i=1:n
    s=0;
    for j=1:N-i+1
        s=s+x(j)*x(j+i-1);
    end
    rr(i)=s/(N-i+1);           % оценка функции корреляции последовательности
end
plot(t,b(1,:),'r',t,rr)
pause
% блоки матрицы точности
D11=[D(1,1) D(1,2) D(1,3) D(1,4)
     D(2,1) D(2,2) D(2,3) D(2,4)
     D(3,1) D(3,2) D(3,3) D(3,4)
     D(4,1) D(4,2) D(4,3) D(4,4)];
D21=[D(5,1) D(5,2) D(5,3) D(5,4)
     D(6,1) D(6,2) D(6,3) D(6,4)];
D22=[D(5,5) D(5,6)
     D(6,5) D(6,6)];
 
C=-inv(D22)*D21                     % матрица связи
BC=inv(D22)                         % корреляционная матрица
[u,v]=eig(BC)
A=u*v^(1/2)*u'                      % оператор окрашивания
 XX=[x(3);x(4);x(5);x(6)]
for i=5:2:N
    xx=randn(2,1)                   % значения белого шума
    y=A*xx+C*XX                     % два следующих значения
    xx(i)=y(1);
    xx(i+1)=y(2);
    XX=[x(i-2);x(i-1);x(i);x(i+1)]        % сдвиг значений 
end
 
 for i=1:n
    s=0;
    for j=1:N-i+1
        s=s+x(j)*x(j+i-1);
    end
    rrr(i)=s/(N-i+1);          % оценка функции корреляции последовательности
end
plot(t,b(1,:),'b',t,rrr)
pause
 
% блоки матрицы точности
D11=[D(1,1) D(1,2) D(1,3)
D(2,1) D(2,2) D(2,3)
D(3,1) D(3,2) D(3,3)];
D21=[D(4,1) D(4,2) D(4,3) 
D(5,1) D(5,2) D(5,3)
D(6,1) D(6,2) D(6,3)];
D22=[D(4,4) D(4,4) D(4,6)
D(5,4) D(5,5) D(5,6)
D(6,4) D(6,5) D(6,6)];
 
CCC=-inv(D22)*D21;           % матрица связи
BCC=inv(D22);                % корреляционная матрица
[u,v]=eig(BCC);
A=u*v^(1/2)*u';             % оператор окрашивания
 
XXX=[x(4);x(5);x(6)]
 
for i=5:3:N
xxx=randn(3,1)              % значения белого шума
y=A*xxx+CCC*XXX 
xxx(i)=y(1);
xxx(i+1)=y(2);
 
XXX=[x(i-1);x(i);x(i+1)]    % сдвиг значений
end
 
 for i=1:n
    s=0;
    for j=1:N-i+1
        s=s+x(j)*x(j+i-1);
    end
    rrrr(i)=s/(N-i+1);      % оценка функции корреляции последовательности
end
plot(t,b(1,:),'g',t,rrrr)