function [c] = savgol(m, nl, nr, ld)
%SAVGOL  Calculate Savitzky-Golay filter coefficients.
%   c = savgol(m,nl,nr,ld)
%   m is the order of the smoothing polynomial (positive integer)
%   nl number of left (past) data points (positive integer)
%   nr number of right (future) data points (positive integer)
%   ld is the order of the derivative desired. ld=0 for smoothing,
%     ld=1 for the filtered first derivate (needs to be divided
%     by the step size). m>=4 is recommended for derivatives.
%
% Example cubic smoothing:
%   f = [1 1.4 1.5 1.4 1.3]'; % example data
%   c = savgol(3, 2, 2, 0);   % smooth with cubic polynomial
%   f_smooth = c'*f           % smoothed value at index = 3
%
% Example calculate first derivative with 2nd order polynomial:
%   f  = [1 1.4 1.5 1.4 1]'; % example data
%   dt = 1.0;                % example data step size
%   c1 = savgol(2, 4, 0, 1); % 1st derivate at rightmost position
%   df_dt = (c1'*f)/dt;      % divide by step size

assert(nl >= 0 && nr >= 0 && ld <= m && nl+nr >= m && ld >= 0);
assert(mod(nl,1)==0 && mod(nr,1)==0 && mod(ld,1)==0 && mod(m,1)==0);

i = -nl:nr;      % rows of Vandermonde matrix
j = 0:m;         % columns of Vandermonde matrix
A = i' .^ j;     % Vandermonde design matrix: A(line, col) = i^j;
K = (A'*A)\A';   % least squares: (A'*A)^(-1)*A'
c = K(ld+1, :)'; % line 1 = smoothing, line 2 = 1st derivative ...

end