function smooth_d=kernel(t,x,h,fun)

%  Kernel Smoothing methods with different kernel functions.
%  Algorithm is from Funcitonal Data Analysis (2nd edition).
%  J.O.Ramsay and B.W.Silverman. Springer series in statistics. 2005.
% 
%  smooth_d=kernel(t,x,h,fun)
%
%  t: sample interval.
% 
%  x: the data for analysis and the data should be a matrix with two columns (rows):
%     one for sample interval and the other for sample value.
% 
%  h: bandwidth parameter.
% 
%  fun: weights function and the default is Quadratic.
%       uni: Uniform.
%       qua: Gasser-Muller weights funciton of Quadratic.
%       fde: Gasser-Muller weights function for the estimation of the first
%            derivative.
%       sde: Gasser-Muller weights function for the estimation of the
%            second derivative.
% 
%  DNP. 2007.12.06


smooth_d=[];
S=[];

if nargin < 3
    error(['Not enough inputs, you should input the data (t and x)' ...
            ' for analysis and choose a bandwidth.']);
    return
elseif nargin >4
    error('Two many inputs!');
    return
elseif nargin==3
    fun='qua';   % Default weight function: Quadratic.  
end

[rown,column]=size(x);
[trown,tcolumn]=size(t);

if ~any([rown column]==1) | ~any([trown tcolumn]==1)
    error('Sample interval and input data should both be a vector!');
    return
elseif length(x)~=length(t)
    error('Number of sample intervals and number of data should be the same!');
    return
elseif rown==1
        x=x';
    end

if all([rown column] <=3 )
    error('Input data should have more than 3 sample points.');
    return
end

if h<=0
    error('Bandwidth h must be a positive value!');
    return
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Kernel Smoothing. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

avg_t=[];

avg_t(1)=t(1);              % average of x with its neighborhood, avg_t(0)=x(1).
avg_t(length(t)+1)=t(end);  % avg_t(end)=x(end).
for i=2:length(t)
    avg_t(i)=(t(i)+t(i-1))/2;
end

for i=1:length(t)
    
    switch fun
        
        case 'uni'
            S=feval('uni',t(i),avg_t,h);
            smooth_d(i)=S*x;
            
        case 'qua'
            S=feval('qua',t(i),avg_t,h);
            smooth_d(i)=S*x;
            
        case 'fde'
            S=feval('fde',t(i),avg_t,h);
            smooth_d(i)=S*x;
            
        case 'sde'
            S=feval('sde',t(i),avg_t,h);
            smooth_d(i)=S*x;
            
    end
    
end
% 
% if rown==1
%     smooth_d=S*x';
% else
%     smooth_d=S*x;
% end



% ******************   Weights functions.  **********************

% ....................

function s=uni(t,avg_t,h)
% Uniform function.

s=0.5*(avg_t(2:end)-avg_t(1:end-1))/h;
s=s.*(abs((avg_t(2:end)-t)/h)<=1);


% ....................

function s=qua(t,avg_t,h)
% Quadratic weight function.

r=(t-avg_t)./h;
s=((3*r(1:end-1)-r(1:end-1).^3)-(3*r(2:end)-r(2:end).^3))./4;
s=s.*(abs((avg_t(2:end)-t)/h)<=1);


% ....................

function s=fde(t,avg_t,h)
% Weigthts for the estimation the first derivative.

r=(t-avg_t)./h;
s=(15*((r(1:end-1).^4-2*r(1:end-1).^2)-(r(2:end).^4-2*r(2:end).^2)))/(16*h);
s=s.*(abs((avg_t(2:end)-t)/h)<=1);


% ....................

function s=sde(t,avg_t,h)
% Weights for the estimation of the second derivative.

r=(t-avg_t)./h;
s=(105*((2*r(1:end-1).^3-r(1:end-1).^5-r(1:end-1))-(2*r(2:end).^3-r(2:end).^5-r(2:end))))/(16*h^2);
s=s.*(abs((avg_t(2:end)-t)/h)<=1);