% the test signal for nonlinear time series with Gaussian statistics will
% be generated
% Author: Sanjay Patil, Oct 31, 2005. Monday. 1012 hours central

% Nonlinear signal model as used in,
% M. Arulampalam, S. Maskell, N. Gordon, and T. Clapp, "A tutorial on
% Particle filters for online Nonlinear/Non-Gaussian Bayesian Tracking,"
% IEEE transactions on signal processing, vol 50, no. 2, February 2002, pp
% 174-188.
 
%
% The set up is described in the section VI of paper, with equations 79-83.
% page 183
% 
% X(k) = Fk(X(k-1), k) + V(k-1)
% Z(k) = Gk(X(k)) + N(k)
% Fk(X(k-1), k) = x(k-1)/2 + 25 * x(k-1)/(1 + x(k-1)^2) + 8 * cos(1.2*k)
% Gk(X(k)) = X(k)^2 / 20
% V(k-1) and N(k) --> zero mean Gaussian random variables with variance
% Q(k-1) and R(k). 
% The paper uses; Q(k-1) = 10, R(k) = 1.

% initial set up
%
close all;
clear all;
clc;

% =========================================================================
% a block to accept inputs from command window (user)
%
% to take the input as a matrix
% Q = input('process noise  variance    : ');
% R = input('measurement noise variance : ');
% x0 = input('initial state value        : ');
% 
% =========================================================================


% declare variables
%
% noise variances
% 
Q = 10;
R = 1;
x0 = 0.2;

len = 40;

% normrnd accepts parameters as;
% mean, standard deviation, 
% enter the value square root of Q in here
%
process_noise = normrnd(0, Q^0.5, 1, len);

% display the parameters of the noise
%
disp(['noise mean     : ',num2str(mean(process_noise))]);
disp(['noise variance : ',num2str(var(process_noise))]);

signal = [];
signal = [x0 signal];
clean_state = x0;
clean_signal = [];
clean_signal = [ x0 clean_signal];
for i = 2:1:len 
    clean_state = clean_state/2 + 25 * clean_state / ( 1 + clean_state^2)...
        + 8 * cos(1.2 * i);
    signal = [signal clean_state];
end;

clean_signal = signal;

%signal = [ signal noise(1:order)];
signal = signal + process_noise;

% display the parameters of the signal
%
disp(['signal mean     : ',num2str(mean(signal))]);
disp(['signal variance : ',num2str(var(signal))]);

% figure of input noise and signal
%
figure;
plot(signal, 'r'); 
hold on; plot(clean_signal, 'k');
% hold on; plot(process_noise); 
title(['Nonlinear Gaussian Time Series - process var = ', num2str(Q)]);
legend('corrupted signal', 'clean signal'); %,'noise');

% to create a additive WGN for observation model
%
measurement = [];
var_measurement_noise = 1.0;
measurement_noise = normrnd(0, var_measurement_noise^0.5, 1, len);
measurement = signal.^2 / 20 + measurement_noise;

% figure of input noise and signal
%
figure;
plot(signal, 'r');hold on; plot(measurement, 'b'); 
hold on; plot(measurement_noise, 'k');
title(['Nonlinear Gaussian Time Series - obsn var = ', num2str(R)]);
legend('corrupted signal', 'measurement', 'measurement noise');

% nlg => nonlinear Gaussian
%
save corrupted_signal_nlg.mat signal
save obsn_nlg.mat measurement
save clean_signal_nlg.mat clean_signal