function [AR, Q, R] = parameterEstimate(frame_data, order)
% function call from kfukfpf.m file - main file line 75
%
% The program to estimate AR parameters for a window of data
%
% input paramters:
% 1. a window of data
% 2. order
%
% output parameters:
% 1. AR coefficients
% 2. process variance
% 3. observation variance
%

% disp('into function parameterEstimate.m line 13');
% % test block for the code to work independently
% %
% load 'obsn.raw';
% frame_data = obsn;
% p = 4;

% The following extracts AR-coefficients for the pseudo-stationary
% analysis-frame
frame = frame_data;
[C Lags] = xcorr(frame, frame);
mid = (length(C) + 1) / 2;

% refer the paper by M. Gabrea
% M. Gabrea, "Adaptive Kalman Filtering-Based Speech Enhancement Algorithm," 
% IEEE conference on audio and speech signal processing, pp I314-I317, 2004.
%
l = 0;
p = order;
C_new = C(mid : mid + 2 * p + l - 1);
R_column = C_new(p + 1 : 2 * p + l);
for i = 1 : p + l
    for j = 1 : p
        R(i,j) = C_new(i + j - 1);
    end;
end;
a = pinv(R) * R_column;   % a(i) is the i'th AR-coefficient

% Computing the observation noise variance using the relation described in
% "Estimation of Noise Variance from Noisy AR signal and it's application
% in speech enhancement", K.K. Paliwal, IEEE transactions on ASSP, vol. 36,
% number 2, Feb., 1998, pp292-294
[C Lags] = xcorr(frame, frame);
mid = (length(C) + 1) / 2;
C_new = C(mid : mid + 2 * p + l - 1);
sum1 = 0;
for i = 1 : p
    sum2 = C_new(i + 1);
    for j = 1 : p
        sum2 = sum2 + a(j) * C_new(abs(i - j) + 1);
    end;
    sum2 = a(i) * sum2;
    sum1 = sum1 + sum2;
    sum1 = sum1 / sum(a.^2);
end;

% measurement noise variance
%
R = sum1;
R = 0.002;


% Defining Process noise variance
% Note, this must be modified to automatically compute the noise
% variance from the data in the analysis frame (using modified Yule-Walker
% equations)
%
% use the AR coefficients to generate back the signal using the first
% 'order' points (assume the first order points to be unaffected)
%
% this is our approach to find the process noise and observation noise
% variance
% 

% % this is causing some trouble - figure it out
% % right now assuming Q to be some multiplier of R
% %
% % 
% noise = normrnd(0, sqrt(R), length(frame), 1);
% x = frame - noise;
% figure;
% plot(frame, 'b');
% hold on;
% plot(x, 'k');
% legend('actual data', 'noise subtracted');
% 
% y = x(1:order);
% 
% for i = order+1:1:length(frame),
%     y(i) = a(1) * y(i-1) + a(2) * y(i-2) + a(3) * y(i-3) + a(4) * y(i-4);
% end;
% 
% % process noise variance
% %
% Q = var(x-y); 
% AR = a;
    
% returning the parameters AR, Q, R.
%
Q = R / 3;
AR = a;

% displaying all the values returned
%
% disp(['The process noise variance     :  ', num2str(Q)]);
% disp(['The measurement noise variance :  ', num2str(R)]);
% disp(['The AR coefficients            :  ']);
% AR

% returning back to respective filter type in main file kfukfpf.m
% 
% disp('into function parameterEstimate.m line 110');
