function [xhat,XHAT,ALPHA]=kallman(U,d,varargin)
%
%[xhat,XHAT,ALPHA]=kallman(x,d,xhat0,lamb,K0,sig)
%Implements unforced kallman filter. Tap weights are the state vectors
%inputs are C(t)
% Copyright (C) 2010 Ikaro Silva
% 
% This library is free software; you can redistribute it and/or modify it under
% the terms of the GNU Library General Public License as published by the Free
% Software Foundation; either version 2 of the License, or (at your option) any
% later version.
% 
% This library is distributed in the hope that it will be useful, but WITHOUT ANY
% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
% PARTICULAR PURPOSE.  See the GNU Library General Public License for more
% details.
% 
% You should have received a copy of the GNU Library General Public License along
% with this library; if not, write to the Free Software Foundation, Inc., 59
% Temple Place - Suite 330, Boston, MA 02111-1307, USA.
% 
% You may contact the author by e-mail (ikaro@ieee.org).

%Initialize all parameters
[N,M]=size(U);
lamb=1;

if (nargin>2 && ~isempty(varargin{1}) )
    xhat=varargin{1};
    xhat=xhat(:);
else
    xhat=rand(M,1)*0.0001;
end
if (nargin>3 && ~isempty(varargin{2}))
    lamb=varargin{2};
end
if (nargin>4 && ~isempty(varargin{3}))
    K=varargin{3};
else
    K=corrmtx(xhat,M-1);
    K=K'*K;
end
if (nargin>5 && ~isempty(varargin{4}))
    sig=varargin{4};
else
    sig=ones(N,M);  %measureement noise (assumed diagonal)!!!!
end


%Do the training
lamb=lamb^(-0.5);
ALPHA=zeros(N,1);
XHAT=zeros(N,M);

for n=1:N
    
    u=U(n,:)';   
    den=diag(1./(u'*K*u + sig(n,:)));
    g= lamb*den*K*u ;
    alpha = d(n) - u'*xhat;
    xhat= lamb *xhat + g*alpha;
    K = (lamb^2)*K - lamb*g*u'*K;

    ALPHA(n)=alpha;
    XHAT(n,:)=xhat';
    
end