% idea to to build a Kalman filter which has the following as the inputs
% 1. Transition Matrix or model, F(n)
% 2. Observation Matrix or model, H(n)
% 3. Initial state vector, x(n)
% 4. Process Noise covariance matrix, Q(n)
% 5. Observation Noise variance, R(n)
% 6. Initial Error Covariance Matrix, P(n)

% let us first start with generating a AR signal
% let us consider the AR signal as used in,
% K. Paliwal, "A noise-compensated long correlation matching method for AR
% spectral estimation of noisy signals," ICASSP 1986, Tokya, 27.7.1-27.1.4
% pp 1369-1372.
% 
% As described in the paper, "The signal is generated from a fourth-order
% AR system (with parameters, a1 = -2.7607, a2 = 3.8106, a3 = -2.6535, a4 =
% 0.9238), driven by zero-mean unit-variance white Gaussian noise process. 
% Generate 100 realizations of the fourth-order AR process and use N=40
% data points during steady-state to compute AR parameters for each of
% these realizations."

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

% =========================================================================
% a block to accept inputs from command window (user)
%
% to take the input as a matrix
% mat1 = input('AR coefficients enter values within rectangular brackets : ');
% 
% =========================================================================


% declare variables
%
% AR filter coefficients
%
a0 = 1;
a1 = -2.7607;
a2 = 3.8106;
a3 = -2.6535;
a4 = 0.9238;

% or following are the values used by Matlab utility
%
% A = [ 1 -2.2137 2.9403 -2.1697 0.9606];

order = 4;
len = 40;
a = [ a0 a1 a2 a3 a4];
var_process_noise = 0.5;
noise = normrnd(0, var_process_noise, len, 1);

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

signal = [];
for i = 1:1:order
    signal(i) = noise(i);
end;

%signal = [ signal noise(1:order)];
for i = 5:1:len
    num = 0;
    for j = 1:1:order
        num = num + a(j+1) * noise(i - j);
    end;
    signal = [signal num];
end;

% 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(noise); hold on; plot(signal, 'r'); legend('noise', 'signal');

% to create a additive WGN for observation model
%
measurement = [];
var_measurement_noise = 2.0;
measurement_noise = normrnd(0, var_measurement_noise, len, 1);
measurement = signal + measurement_noise';

% figure of input noise and signal
%
figure;
plot(signal, 'r');hold on; plot(measurement, 'b'); 
legend('signal', 'measurement');

save signal.mat signal
save obsn.mat measurement

