% function kfukfpf (filterType, frame_size, window_size, fileName, samplingFrequency, order)
% main block of the filter code 
% Author : Sanjay Patil, Nov 02, 2005. Wednesday. 1536 hours central
%
%
% function to run the filtering algorithm given the following parameters
% input paramters:
% 1.filterType - 
%          1 -> Kalman filter, 
%          2 -> unscented Kalman filter, 
%          3 -> particle filter
% 2. frame_size - frame size for which the signal is cleaned (ms)
% 3. window_size - window size for which the model parameters are estimated
% 4. fileName - the binary raw file to be processed
% 5. samplingFrequency - sampling frequency as KHz
% 6. order - AR order - linear model considered
% frame size and window size are milliseconds (ms)

% check for the filter type
%
% if filterType == 'KFI'
%     disp('Kalman filter implementation');
%     disp(['the character string is ', filterType]);
% elseif filterType == 'UKF'
%     disp('Unscented Kalman filter implementation');
%     disp(['the character string is ', filterType]);
% elseif filterType == 'PFI'
%     disp('Particle filter implementation');
%     disp(['the character string is ', filterType]);
% end;

% clear all;
close all;
clc;

% initial values
%
filterType = 'KFI';
frame_size = 70;
window_size = 70;
fileName = 'obsn_ar.raw';
samplingFrequency = 8;
order = 4;

% right now assume that the file is single frame only
%
fid = fopen(fileName, 'r');
x = fread(fid, 'int16');
fclose(fid);

% x = y(1650:5600);
% figure;
% plot(x); title('input signal');
% range = max(x) - min(x);

x = x / 2^15; % range;

figure; plot(x); title('normalized input signal line 58 kfukfpf.m');

frame_len = samplingFrequency * frame_size;
window_len = samplingFrequency * window_size;

% The code considers left aligned frame for a window analysis 
% refer to http://www.cavs.msstate.edu/hse/ies/projects/speech/software
% /documentation/class/algo/Window/index.html
% This analysis is significant for overlapping of the windows
% 
if frame_len < window_len
    len = length(x) - window_len;
    num_frames = floor(len/frame_len);
elseif frame_len == window_len
    len = length(x) - window_len;
    num_frames = floor(len/frame_len);
elseif frame_len > window_en
    len = length(x) - frame_len;
    num_frames = floor(len/window_len);
end;

parameter_estimate_block_len = window_len;
block_jump = frame_len;
clean_signal = [];
for frame_index = 0:1:num_frames
    

    % % this will not work right now, as frame_data is not declared
    % % to display the input signal here
    % figure (1);
    % plot(frame_data, 'r');
    % legend('corrupted speech signal - one frame only');

    % call the proper function now
    % depending on the filterType

    if filterType == 'KFI'
        frame_data = x(frame_index * block_jump + 1 : frame_index * block_jump + frame_len);

        % figure; plot(frame_data);
        % title([num2str(frame_index+1), 'th frame of data']);

        % estimating the AR coefficients and process and observation noise
        %
        window_analysis_frame = x(frame_index * block_jump + 1 : frame_index * block_jump + parameter_estimate_block_len);
        [AR, Q, R] = parameterEstimate(window_analysis_frame, order);
        if frame_index == 0
            % assign the first order values as clean speech signal itself
            %
            initial_state = frame_data(1:order);
            clean_signal = kalmanFilter(AR, Q, R, frame_data, initial_state, frame_index);
        else
            % now the initial state will be the clean data
            %
            disp(['the frame code is working on : ', num2str(frame_index+1)]);
            initial_state = clean_signal(end-order+1:end);
            filtered_signal = kalmanFilter(AR, Q, R, frame_data, initial_state, frame_index);
            clean_signal = [clean_signal filtered_signal];
        end;
    else
        disp('code not implemented yet!');
    end;
 
    figure;
    plot(x(1:length(clean_signal))); hold on; plot(clean_signal, 'r');
    pause;
end;

figure;
x = x;
clean_signal = clean_signal;
plot(x); hold on; plot(clean_signal, 'r');
legend('input signal', 'output signal');
title('unnormalized signal kfukfpf.m line 131');

figure;
range_x = max(x) - min(x);
x = x / range_x;
range_clean = max(clean_signal) - min(clean_signal);
clean_signal = clean_signal / range_clean;
plot(x); hold on; plot(clean_signal, 'r');
legend('input signal', 'output signal');

% just a debugging idea here
% will be deleted afterwards
% we know the name of the clean signal file
% we will use that file to find the MSE and put it here
%
fid1 = fopen('clean_ar.raw', 'r');
data = fread(fid1, 'int16');
data = data / 2^15;
clean_signal = clean_signal * range_clean;
MSE = mean(square(clean_signal-data))