% function signal_generator_ar(AR_model, Q, R, initial_value)
% this is a module that will accept the AR coefficients, order, 
% noise variance (process), and noise variance (observation) 
% and store the clean and the observation signal 
% 
% Author: Sanjay Patil Nov 02, 2005. Wednesday. 1018 hours central
% last modified : Nov 03, 2005. Thursday. 1748 hours central
%

% this will the command line argument for the signal generation afterwards
%
% function signal_generator_ar(AR_model, Q, R, initial_value)

% if you are having a function all, do not free memory allocated to the 
% variables
%
% initial set-up
%
% clear all;
close all;
clc;

% details of the program
%
disp('The program to generate test signal for the AR model');
disp('Input Parameters (four parameters) :');
disp('1. AR coefficients  2. Process noise variance ');
disp('3. Observation noise variance, 4. initial values : ');
disp('Output parameters (two files) :');
disp( '1. clean file (clean.raw)');
disp('2. obsn file (obsn.raw)');

disp(' ');

% % uncomment afterwards
% %
% disp('displaying all the values :');
% 
% disp(['the AR Coefficients            : ', num2str(AR_model)]);
% disp(['the Process noise variance     : ', num2str(Q)]);
% disp(['the Observation noise variance : ', num2str(R)]);
% disp(['the initial value              : ', num2str(initial_value)]);
% 
% 
% disp('displayed all the values :');

% % accepting values from command window
% %
% AR_model = input('Enter the array of AR coefficients [row matrix] : ');
% Q = input('Enter the process noise variance : ');
% R = input('Enter the observation noise variance : ');
% initial_value = input('Enter the initial values [row matrix] : ');

% example values : 
Q = 0.002;
R = 0.005;
AR_model = [0.03 0.06 0.70 0.05];
initial_value = [0.1 0.2 0.5 0.8];

x(1) = initial_value(1);    z(1)= x(1) + normrnd(0, sqrt(R));
x(2) = initial_value(2);    z(2)= x(2) + normrnd(0, sqrt(R));
x(3) = initial_value(3);    z(3)= x(3) + normrnd(0, sqrt(R));
x(4) = initial_value(4);    z(4)= x(4) + normrnd(0, sqrt(R));
 
% The state space model is a 4'th order AR process 
% with an AWGN (0 mean, 0.02 variance) driving noise. 
% The measurement model has signal added with obsn noise variance.
for i = 5 : 560
    x(i)=AR_model(1) * x(i-1) + AR_model(2) * x(i-2) +...
        AR_model(3) * x(i-3) + AR_model(4) * x(i-4) +...
        normrnd(0, sqrt(Q));
    z(i)= x(i) + normrnd(0, sqrt(R)); %Adjust the noise variance here
end;

% save ('clean.raw', 'x', '-ascii', '-double');
% save ('obsn.raw', 'z', '-ascii', '-double');

figure;
plot(x, 'b');
hold on;
plot(z, 'm');
title(['Time series : AR model, process noise N(0,{\sigma} = ',...
    num2str(Q), ') and measurement noise N(0,{\sigma} = ', num2str(R), ']']);
legend('clean signal', 'observation signal');

x = x * 2^15;
z = z * 2^15;
fid1 = fopen('clean_ar.raw', 'w');
fwrite(fid1, x, 'int16');
fclose(fid1);

fid2 = fopen('obsn_ar.raw', 'w');
fwrite(fid2, z, 'int16');
fclose(fid2);
 
% testing back the results here itself
%
x = x / 2^15;
z = z / 2^15;
r = xcorr(z);
r1 = xcorr(x);
r(1:length(x)-1) = []; 
r1(1:length(z)-1) = [];
ar_coef_x = levinson(r,4)
ar_coef_z = levinson(r1,4)
ar_coef_x = -ar_coef_x(2:end);
ar_coef_z = -ar_coef_z(2:end);

% we will try to reconstruct the signal back here again
% using the same set up
% let us assume the same Q and R values and go ahead
% 
% reconstructed_x = [];
reconstructed_x = x(1:4);
% reconstructed_z = [];
reconstructed_z = z(1:4);

for i = 5:length(x)
    % the corrupted or observed signal
    %
    reconstructed_z(i) = ar_coef_z(1) * reconstructed_z(i-1) + ...
       ar_coef_z(2) * reconstructed_z(i-2) + ...
       ar_coef_z(3) * reconstructed_z(i-3) + ... 
       ar_coef_z(4) * reconstructed_z(i-4) + normrnd(0, sqrt(R));

   % the original (uncorrupted) signal
   %
   reconstructed_x(i) = ar_coef_x(1) * reconstructed_x(i-1) + ...
       ar_coef_x(2) * reconstructed_x(i-2) + ...
       ar_coef_x(3) * reconstructed_x(i-3) + ... 
       ar_coef_x(4) * reconstructed_x(i-4) + normrnd(0, sqrt(Q));
end;

figure;
plot(reconstructed_x, 'b');
hold on;
plot(reconstructed_z, 'm');
title(['Reconstructed Time Series : AR model, process noise N(0,{\sigma} = ',...
    num2str(Q), ') and measurement noise N(0,{\sigma} = ', num2str(R), ']']);
legend('Resconstructed clean signal', 'Reconstructed observation signal');

   
   