% Kalman filtering of speech signals
% This code implements a Kalman filter on the acoustic samples of speech.
% Authors: |- Saurabh Prasad
%          |- Sanjay Patil

clear;
[y fs] = wavread('bb_1315a.wav');
p = 8;  % Define the AR order
l = 0;
frame_length = 80;
x = y(14000:14000+frame_length);

x_corrupted = x + normrnd(0,0.02,length(x),1); % Corrupt the acoustic space 
                                              % with AWGN (0 mean, 0.2 var.)

% The following extracts AR-coefficients for the pseudo-stationary
% analysis-frame
[C Lags] = xcorr(x_corrupted,x_corrupted);
mid = (length(C)+1)/2;
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

% The a(i)'s are then used to generate the state-transition matrix 
% representing state evolution in that frame
F_ = eye(p-1,p-1);
F = zeros(p,p);
F(:,1) = zeros(p,1);
F(p,:) = a';
% F is now the state-transitCleanion matrix for the Kalman Filter
F(1:p-1,2:p) = eye(p-1,p-1);

% Defining G (Process noise gain) and H (observation matrix)
G = zeros(p,1);
G(p,1) = 1;
H = G';

% Initialize the initial state-vector prediction to the first 'p' data
% samples
X(:,1) = y(1:p);

% 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(x_corrupted,x_corrupted);
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;
sigma_v = sum1;

% Defining Process noise variance
sigma_u = .15;

P_old = eye(p,p);
Q = sigma_u*eye(p,p);
R = sigma_v;

s(1:p) = x_corrupted(1:p); % not filtering the first p data points in the analysis frame
x_state = x_corrupted(1:p);
for i = p+1 : frame_length
   z = x_corrupted(i);
%   x_state = s(i-p:i-1)'; % The state vector passed to the Kalman filter must be px1
   [s(i) temp] = kalman(x_state,z,Q,R,F,G,H);   
   x_state = temp;
end;

hold on; plot(x_corrupted,'r');
hold on; plot(s);
hold on; plot(x,'g');

     
 
    