hnstm=sgsdf_gram_poly(n,s,t,m)

function hnstm=sgsdf_gram_poly(n,s,t,m)
% Savitzky-Golay smoothing and differentiation filter 
% (i.e.,polynomial smoothing and differentiation filte, 
%  or least-squares smoothing and differentiation filter) 
% calculated with the recursive properties of Gram polynomials. 
% This reccursive procedure can avoid the inaccuracy 
% due to singular or badly scaled matrix when using a matrix approach.
%
% hnstm=sgsdf_gram_poly(n,s,t,m)
% n:      polynomial degree
% s:      derivative(differentiation) order (0=smoothing)
% t:      evaluation point (commonly,t=0,i.e.,smoothing or differentiation
%               at the central point of the 2*m+1 points)
% 2*m+1:  data point number
% hnstm:  convoluction weights  
%
% Examples:
% Table I in ref.[2] 
% analytic solution
% for m=sym(2:3)
%     n=sym(2);
%     for s=sym(0:1)
%         disp(sprintf('Savitzky-Golay Smoothing and Differentiation Filter(%d points,%d-th polynomial,%d-derivatie)',double(2*m+1),double(n),double(s)));
%         for t=-m:m
%             disp(sgsdf_gram_poly(n,s,t,m));
%         end
%     end
% end
%
% Table II in ref.[2] 
% numerical solution (much faster than analytic solution) 
% for m=2:3
%     n=3;
%     for s=0:1        
%         disp(sprintf('Savitzky-Golay Smoothing and Differentiation Filter (%d points,%d-th polynomial,%d-derivatie)',2*m+1,n,s));
%         for t=-m:m            
%             disp(sgsdf_gram_poly(n,s,t,m));
%         end
%     end    
% end
%
% Table III in ref.[2] 
% numerical solution
% for m=2:10
%     n=2;
%     s=0;
%     t=-m;
%     disp(sgsdf_gram_poly(n,s,t,m));
% end
%
% Table IV in ref.[2] 
% numerical solution 
% for m=2:10
%     n=2;
%     s=1;
%     t=-m;
%     disp(sgsdf_gram_poly(n,s,t,m));
% end
%
%Author:
%       Jianwen Luo,Ph.D Candidate
%       Email:[e-mail], [e-mail]
%       Date: July 25,2004
%       Department of Biomedical Engineering
%       Tsinghua University, Beijing 100084, P. R. China 
%
% References:
% [1]A. Savitzky and M. J. E. Golay, "Smoothing and Differentiation of Data by Simplified Least Squares Procedures," 
%   Analytical Chemistry, vol. 36, pp. 1627-1639, 1964.
% [2]P. A. Gorry, "General Least-Squares Smoothing and Differentiation by the Convolution (Savitzky-Golay) Method," 
%   Analytical Chemistry, vol. 62, pp. 570-573, 1990.
%
%See also sgolay (Matlab function)
%         sgsdf  (http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=4038&objectType=file)   

hnstm=[];
% i:      convolution weight of the i-th point (-m<=i<=m)
for i=-m:m
    hnstm=[hnstm hinstm(i,n,s,t,m)];
end

function value=hinstm(i,n,s,t,m)
value=0;
for k=0:n
    value=value+(2*k+1)*gff(2*m,k)/gff(2*m+k+1,k+1)*Pkmi(k,m,i)*Pkmsi(k,m,s,t);%(see ref.[2])
end
 
function value=gff(a,b)
%generalized factorial function (see ref.[2])
if b==0
    value=1;
else
    value=prod(a-(0:b-1));
end

function value=Pkmi(k,m,i)
%Gram polynomials (equ.(10) in ref.[2])
if k==-1
    value=0;
elseif k==0
    value=1;
else
    value=2*(2*k-1)/k/(2*m-k+1)*i*Pkmi(k-1,m,i)-(k-1)*(2*m+k)/k/(2*m-k+1)*Pkmi(k-2,m,i);
end

function value=Pkmsi(k,m,s,i)
%s-th derivative of Gram polynomials (equ.(11) in ref.[2])
if s==0
    value=Pkmi(k,m,i);
else    
    if k==-1
        value=0;
    elseif k==0
        value=0;
    else
        value=2*(2*k-1)/k/(2*m-k+1)*(i*Pkmsi(k-1,m,s,i)+s*Pkmsi(k-1,m,s-1,i))-...
            (k-1)*(2*m+k)/k/(2*m-k+1)*Pkmsi(k-2,m,s,i);       
    end   
end