% 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.001;
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_ar1.raw', 'w');
fwrite(fid2, z, 'int16');
fclose(fid2);
 
% testing back the results here itself
%
x = x / 2^15;
z = z / 2^15;
r_z = xcorr(z);
r_x = xcorr(x);
r_x(1:length(x)-1) = []; 
r_z(1:length(z)-1) = [];
disp('clean AR signal AR coefficients : ');
ar_coef_x = levinson(r_x,4)
disp('corrupted AR signal AR coefficients : ');
ar_coef_z = levinson(r_z,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');

% the observation here is:
% one example run, that was made
% clean AR signal AR coefficients : ar_coef_x =
%     1.0000   -0.1656   -0.1385   -0.6550    0.0617
% corrupted AR signal AR coefficients : ar_coef_z =
%     1.0000   -0.1675   -0.1556   -0.4975   -0.0062

% let me try to find the ar coefficients from the method as described in M.
% Gabrea, "Adaptive Kalman filtering technique for speech enhancement,"
% ICASSP. [please the name of the paper].
% page 523.
%
% [LPC] = [ R_obsn p by p+l]' * [R_obsn 1 by p+l] equation 15 of the paper
%
r_z_icassp = xcorr(z);
p = length(AR_model);
off_center_index = find(r_z_icassp == max(r_z_icassp)) + 1;

r_z_icassp(find(r_z_icassp == max(r_z_icassp)));
r_z_icassp(find(r_z_icassp == max(r_z_icassp))+1);

l = 1;
ar_len = 2 * p + l;
r_z_icassp = r_z_icassp(off_center_index:end);
r_z_icassp = r_z_icassp(1:ar_len);

mat_z = [];
for i = 1:p+l
    for j = 1:p
        mat_z(i,j) = r_z_icassp(i+j-1);
    end;
end;


vec_z = r_z_icassp(p+1:ar_len)';

Lpc_z = pinv(mat_z) * vec_z;

Lpc_z

r_x_icassp = xcorr(x);
off_center_index = find(r_x_icassp == max(r_x_icassp)) + 1;

r_x_icassp(find(r_x_icassp == max(r_x_icassp)));
r_x_icassp(find(r_x_icassp == max(r_x_icassp))+1);

l = 1;
ar_len = 2 * p + l;
r_x_icassp = r_x_icassp(off_center_index:end);
r_x_icassp = r_x_icassp(1:ar_len);

mat_x = [];
for i = 1:p+l
    for j = 1:p
        mat_x(i,j) = r_x_icassp(i+j-1);
    end;
end;


vec_x = r_x_icassp(p+1:ar_len)';

Lpc_x = pinv(mat_x) * vec_x;

% Lpc_x  


% now to compute the noise variance based on K. Paliwal, "Estimation of
% noise variance from the noisy AR signal and its application in speech
% enhancement," IEEE, 1988, pp 292-294.

%  step 1. generate autocorrelation coefficients Ry(i), i = 0, 1... p+q,
%  where q > p
%  step 2. compute the AR coefficients using the above method
%  step 3. using the above AR coefficients and the equation (6) from the
%  referred paper.
%
% step 1. generate the autocorrelation coefficients Ry(i), i = 0, 1,... p+q
% where q > p
%
% figure;
% plot(x, 'b');
% hold on; 
% plot(z, 'k');
% legend('clean signal','corrupted signal');


r_z_ieee = xcorr(z);
off_center_index = find(r_z_ieee == max(r_z_ieee));

q = p + 1;
ar_len = p + q;
r_z_ieee = r_z_ieee(off_center_index : off_center_index + ar_len);
r_z_ieee = r_z_ieee / (max(r_z_ieee) - min(r_z_ieee));

% step 2. compute the AR coefficients using the above method
%
mat_z_ieee = [];
for i = 1:q
    for j = 1:p
        mat_z_ieee(i,j) = r_z_ieee(i+j);
    end;
end;

vec_z_ieee = r_z_ieee(p+2:ar_len + 1)';
Lpc_z_ieee = pinv(mat_z_ieee) * vec_z_ieee;
Lpc_z_ieee

% step 3. find the process noise variance using equation(6)
%

denominator = sum(Lpc_z_ieee.^2);
noise_var = 0;
for i = 1:1:p
    sum_inner = 0;
    for k = 1:1:p
        sum_inner = sum_inner + Lpc_z_ieee(i) * r_z_ieee(abs(i-k)+1);
    end;
    sum_inner = sum_inner + r_z_ieee(i+1);
    noise_var = noise_var + sum_inner * Lpc_z_ieee(i);
end;

noise_var = noise_var / denominator;
noise_var


% =======================================================================
%                               RESULTS
% =======================================================================

% LPC computed based on the above referred paper for the different runs 
% this does indicate the comment made in the paper, that "even if this
% approach may sometimes lead to unsatisfactory performances," still the
% paper proposes to estimate the AR parameters from the modified
% Yule-Walker equations
% For example, refer to the results gotten for run # 6.
%


% % run # 1
% Lpc_z =     0.3147     0.7169     0.1655    -0.4650
% Lpc_x =    -0.1157     0.7029     0.0340     0.2529

% % run # 2
% Lpc_z =    -0.3118     0.6256    -0.0398     0.6659
% Lpc_x =     0.3114     0.6698     0.0960    -0.2911

% % run # 3
% Lpc_z =    -0.1958     0.2474     0.0951     0.7914
% Lpc_x =    -0.2275     0.6459    -0.0051     0.5011

% % run # 4
% Lpc_z =    -0.1721     0.8369    -0.1736     0.4516
% Lpc_x =    -0.3186     0.7453    -0.0353     0.5394

% % run # 5
% Lpc_z =    -0.3636     0.6616     0.1914     0.4559
% Lpc_x =    -0.5656     0.5006     0.0532     0.9936

% % run # 6
% Lpc_z =     3.7155    0.8112     -0.9172    -2.8264
% Lpc_x =     0.2555     0.7508     0.0122    -0.2008

% % run # 7
% Lpc_z =     0.2028     0.7985    -0.1650     0.0448
% Lpc_x =    -0.8462     0.6958    -0.0518     1.2217

% % run # 8
% Lpc_z =     -0.1830     0.8750     0.2390    -0.0286
% Lpc_x =     -0.0197     0.7529     0.0841     0.0708
