Приложение 1

Текст m-файла к примеру выполнения работы №2

clear all; close all;

%Исходные данные

s0 = 1;

s1 = 5;

s2 = 2;

r = 50;

T = 5;

m = 50;

dt = T/m; %период дискретности

%Коэффициенты полинома

x = normrnd(0, [s0; s1; s2]);

x0 = x(1); x1 = x(2); x2 = x(3);

%Моделирование измерений

R = r^2*eye(m);

t = 0:dt:T-dt;

H = [ones(m,1) t' t'.^2];

v = normrnd(0,r,m,1);

y = H*x + v;

%Расчет оценок МНК

xmnk = (H'*H)^-1*H'*y

Pmnk = (H'*H)^-1*H'*R*H*(H'*H)^-1;

eps = x - xmnk

sigma = sqrt(diag(Pmnk))

ymnk = H*xmnk;

%Расчет линейных байесовских оценок

P0 = diag([s0^2 s1^2 s2^2]);

Pb = (P0^-1 + H'*R^-1*H)^-1;

xb = Pb*H'*R^-1*y

eps = x - xb

sigma = sqrt(diag(Pb))

yb = H*xb;

%Построение графиков

figure;

plot(t,y);

xlabel('t,c');

ylabel('y(t)');

title('Измерения');

figure;

plot(t,H*x,'b--',t,ymnk,'r',t,yb,'k');

xlabel('t,c');

ylabel('y~(t)');

title('Оценки тренда');

legend('истинный тренд','оценка МНК','лин. байесовская оценка');

Приложение 2

Текст m-файла к примеру выполнения работы №3

clear; close all;

N = 200000; %число отсчетов случайного процесса

M = N/200; %число отсчетов корреляционной ф-ии

%Параметры случайного процесса

sig = 1.5;

alf = 0.5;

bet = 2;

dt = 0.01; %период дискретности

%Формирование непрерывной модели динамики

F = [0 1; -(alf^2 + bet^2) -2*alf];

G = [0; sig*sqrt(2*alf)];

H = [sqrt(alf^2 + bet^2) 1];

%Дискретизация непрерывной модели

Fi = eye(2) + F*dt;

Gam = G*dt;

%Формирование начального состояния фильтра

P = [sig^2/2/(alf^2 + bet^2) 0; 0 sig^2/2];

x = sqrt(P)*normrnd(0,1,2,1);

%Формирующий фильтр

y(N) = 0;

w = normrnd(0,1/sqrt(dt),1,N);

for k = 1:N

x = Fi*x + Gam*w(k);

y(k) = H*x;

end;

%Обработка результатов моделирования

t = (-M:M)*dt;

%теоретическая корреляционная функция

r2 = sig^2*exp(-alf*abs(t)).*cos(bet*t);

%выборочная корреляционная функция

R2 = xcov(y,'unbiased');

m1 = 1:N/100;

figure;

plot(m1*dt,y(m1));

xlabel('t, s');

title('Stochastic process');

figure;

plot(t,R2(N-M:N+M));

hold on;

plot(t,r2,'r--');

xlabel('\tau, s');

ylabel('R(\tau)');

title('Correlation function');

Приложение 3

Текст m-файла к примеру выполнения работы №4

clear; close all;

N = 1000; %число отсчетов случайного процесса

%Параметры случайного процесса

sig = 1.5;

alf = 0.5;

bet = 2;

sigv = 0.5; %ско измерений

dt = 0.01; %период дискретности

%Формирование непрерывной модели динамики

F = [0 1; -(alf^2 + bet^2) -2*alf];

G = [0; sig*sqrt(2*alf)];

H = [sqrt(alf^2 + bet^2) 1];

H1 = H;

%Дискретизация непрерывной модели

Fi = eye(2) + F*dt;

Gam = G*dt;

Q = 1/dt;

%Формирование начального состояния фильтра

P = [sig^2/2/(alf^2 + bet^2) 0; 0 sig^2/2];

x = sqrt(P)*normrnd(0,1,2,1);

xo = [0; 0]; %начальная оценка

%Инициализация массивов

%для записи результатов моделирования

y(N) = 0;

Xo(2,N) = 0;

sko(2,N) = 0;

X(2,N) = 0;

%Формирующий фильтр

w = normrnd(0,1/sqrt(dt),1,N); %порождающий шум

v = normrnd(0,sigv,N,1); %шум измерений

for k = 1:N

x = Fi*x + Gam*w(k);

y(k) = H*x + v(k);

X(:,k) = x;

end;

%обработка измерений

R = sigv^2;

Pe = P;

for k=1:N

if (k==2/dt || ((k>=4/dt) && (k<=7/dt)))

H = H1; %Получение измерений

else H = H1*0;

end;

xp = Fi*xo; %Прогноз

S = Fi*Pe*Fi' + Gam*Q*Gam'; %Ковариация прогноза

K = S*H'*(H*S*H' + R)^-1; %К.у. фильтра Калмана

xo = xp + K*(y(k) - H*xp); %Оценка

Pe = (eye(2) - K*H)*S; %Ковариация ошибок

Xo(:,k) = xo;

sko(:,k) = sqrt(diag(Pe));

end;

%Построение графиков ошибок оценивания

t = (0:N-1)*dt;

figure;

plot(t,X(1,:) - Xo(1,:));

hold on;

plot(t,3*sko(1,:),'r--');plot(t,-3*sko(1,:),'r--');

xlabel('t, s');

ylabel('\epsilon_1');

title('Error of x_1 estimation');

figure;

plot(t,X(2,:) - Xo(2,:));

hold on;

plot(t,3*sko(2,:),'r--');plot(t,-3*sko(2,:),'r--');

xlabel('t, s');

ylabel('\epsilon_2');

title('Error of x_2 estimation');

39