compare_error_variance_monte_carlo.m

% Error variance analysis for 100 simulation runs (Monte Carlo)
%
% An example for Rapid STM32 Blockset
%
% Application Note AN1 Developing a Tilt Sensor System Using Rapid STM32 Blockset
%
% Copyright 2009 Krisada Sangpetchsong
%
% Visit  www.rapidstm32.com for further information

clc % clear screen

start = 0; % sec
stop = 10; % sec
final_run = 100; % Number runs to simulate

% Generate random seeds to drive the Band-Limited White Noise Block in
% the Accelerometer Sensor Subsystem
% the seed is picked from a uniform random variable between begin_bound and
% end_bound
RandStream.setDefaultStream(RandStream('mt19937ar','seed',sum(100*clock))); % set random seed based on clock
begin_bound = 0;
end_bound = 10000;
seed_vec = floor(begin_bound + (end_bound - begin_bound).*rand(2,final_run));

for k = 1:final_run
    seed1 = seed_vec(1,k);
    seed2 = seed_vec(2,k);
    
    disp(['Run' num2str(k)])
    sim('tilt_system',[start stop]); % This line of code starts and stops Simulink simulation automatically.
    
    if (k == 1)
        % First run
        % initialise record buffer
        true_roll_error_deg_rec = repmat(error_variance.signals(1).values(:,1),1,final_run);
        theoretical_roll_error_deg_rec = repmat(error_variance.signals(1).values(:,2),1,final_run);
        true_pitch_error_deg_rec = repmat(error_variance.signals(2).values(:,1),1,final_run);
        theoretical_pitch_error_deg_rec = repmat(error_variance.signals(2).values(:,2),1,final_run);
    else
        % Other runs
        % Import data from saved work space data
        true_roll_error_deg_rec(:,k) = error_variance.signals(1).values(:,1);
        theoretical_roll_error_deg_rec(:,k) = error_variance.signals(1).values(:,2);
        true_pitch_error_deg_rec(:,k) = error_variance.signals(2).values(:,1);
        theoretical_pitch_error_deg_rec(:,k) = error_variance.signals(2).values(:,2);
    end
end

time = error_variance.time;
mean_true_roll_error_deg = mean(true_roll_error_deg_rec,2);
mean_theoretical_roll_error_deg = mean(theoretical_roll_error_deg_rec,2);
mean_true_pitch_error_deg = mean(true_pitch_error_deg_rec,2);
mean_theoretical_pitch_error_deg = mean(theoretical_pitch_error_deg_rec,2);

std_true_roll_error_deg = std(true_roll_error_deg_rec,0,2);
std_true_pitch_error_deg = std(true_pitch_error_deg_rec,0,2);

% Plot results
subplot(2,1,1)
plot(time,mean_true_roll_error_deg,'r-',time,mean_theoretical_roll_error_deg,'b-',time,...
    std_true_roll_error_deg,'g-')
ylabel('Roll Errors (deg)'), xlabel('time (sec)')
legend('Mean errors','\surd(P_{11})','Error std')

subplot(2,1,2)
plot(time,mean_true_pitch_error_deg,'r-',time,mean_theoretical_pitch_error_deg,'b-',time,...
    std_true_pitch_error_deg,'g-')
ylabel('Pitch Errors (deg)'), xlabel('time (sec)')
legend('Mean errors','\surd(P_{22})','Error std')