function varargout = sg(varargin)
% 
%
%   Filter creating:
%   CM = sg(w, p, m)  
%
%   Filter applying:
%   Y_hat = sg(CM, Y) 
%   Y_hat = sg(CM, Y, dim) 
%
%   Filter creating & applying:
%   [Y_hat, CM] = sg(w, p, m, Y)
%   [Y_hat, CM] = sg(w, p, m, Y, dim)
%
%
% w     window size in points; w must be:  (odd & 3 <= w < d) or (w == d); 
%       (d is the size of the dimension of the input matrix Y along which it should be filtered,
%       if Y is not supplied, it is assumed that d == inf)
%
% p     polynom order p=0,1,2,...
%
% m     derivative order m=0,1,2...
%
% Y     matrix which should be smoothed/differinate; 
%       by default (i.e. if the parameter dim is not supplied) 
%       filtering will be applied to the first non-singleton dimension; 
%
% dim   selects the dimension along which Y should be filtered
%
%
% CM    coefficient matrix (filter) 
% Y_hat result matrix
%

% Copyright: S.Verzakov, serguei@ph.tn.tudelft.nl
% Faculty of Applied Sciences, Delft University of Technology
% P.O. Box 5046, 2600 GA Delft, The Netherlands

% $Id: sg.m,v 1.9 2007/03/15 16:10:24 serguei Exp $

	if nargout > 0
		varargout(1:nargout) = {[]};
	end
	
	Y_hat  = [];
	CM     = [];

	w   = [];
	p   = [];
	m   = [];
	Y   = [];
	dim = [];

	d   = [];

	create_filter = 0;
	apply_filter  = 0;

	switch length(varargin)
	 case 2
	  %Y_hat = sg(CM, Y) 
	  apply_filter  = 1;
	  CM  = varargin{1};
	  Y   = varargin{2};

	 case 3
	  %CM = sg(w, p, m)  
	  if all(size(varargin{1}) == 1)
		  create_filter  = 1;
		  w  = varargin{1};
		  p  = varargin{2};
		  m  = varargin{3};

		  %Y_hat = sg(CM, Y, dim) 
	  else
		  apply_filter  = 1;
		  CM  = varargin{1};
		  Y   = varargin{2};
		  dim = varargin{3};
	  end  

	 case 4
	  %[Y_hat, CM] = sg(w, p, m, Y)
	  create_filter = 1;
	  apply_filter  = 1;
	  w  = varargin{1};
	  p  = varargin{2};
	  m  = varargin{3};
	  Y  = varargin{4};

	 case 5
	  %[Y_hat, CM] = sg(w, p, m, Y, dim)
	  create_filter = 1;
	  apply_filter  = 1;
	  w   = varargin{1};
	  p   = varargin{2};
	  m   = varargin{3};
	  Y   = varargin{4};
	  dim = varargin{5};

	 otherwise
	  error('Invalid number of input parameters');
	end


	if apply_filter
		if ~isempty(dim) 
			if floor(dim) ~= dim | dim < 1 | dim > length(size(Y)) 
				error('dim must be integer, >=1, <=length(size(Y))');
			end  
		else 
			dim = size(Y);
			dim = find(dim > 1);

			if isempty(dim)
				error('Filter cannot be applied to the scalar Y');
			end

			dim = dim(1);
		end

		d = size(Y,dim);

		if d < 3
			error('size(Y,dim) must be >= 3');
		end
	end


	if create_filter
		if ~isempty(d) 
			if floor(w) ~= w | mod(w,2) == 0 | w < 3 
				error('w must be odd integer >=3');
			end
		else
			if (floor(w) ~= w | mod(w,2) == 0 | w < 3 | w > d) & w ~= d
				error('w must be (odd integer >=3, <=size(Y,dim)) or ==size(Y,dim)');
			end
		end  
		
		M = floor(w/2);

		if floor(p) ~= p | p < 0
			error('p must be 0,1,2...');
		end

		if floor(m) ~= m | m < 0
			error('m must be 0,1,2...');
		end

		if m > p
			Y_hat = zeros(size(Y));
			CM = zeros(3);

			WarnMsg = 'Derivative order is greater than order of approximating polinomial';
			if exist('prwarning')
				prwarning(1, WarnMsg);
			else
				warning(WarnMsg);
			end
			
		else
			X = repmat([(-M):(w-M-1)]',[1 p]);
			X = [fliplr(cumprod(X,2)) ones(w,1)];

			CM = (X'*X)\X';
			%CM = X\eye(w);
			CM = CM(1:p-m+1,:);

			if m == 0
				CM = X*CM;
			else
				factors = zeros(1,p-m+1);
				for i=m:p
					factors(p-i+1) = prod([1 (i-m+1):i]);
				end
				CM = [X(:,(m+1):(p+1))].*repmat(factors,[w 1])*CM;
			end
		end
	end


	if apply_filter & isempty(Y_hat)
		cms = size(CM);

		if cms(1) ~= cms(2) 
			error('size(CM,1) must be equal size(CM,2)');
		elseif (mod(cms(2),2) == 0 | cms(2) < 3 | cms(2) > d) & cms(2) ~= d
			error('size(CM,2) must be (odd, >=3, <=size(Y,dim)) or == size(Y,dim)');
		end

		w = cms(2);
		M = floor(w/2);

		ys = size(Y);
		order = [1:length(ys)];
		order(1) = dim;
		order(dim) = 1;
		Y = permute(Y,order);
		ysn = size(Y);

		CM_b = CM(1:M,:);
		CM_m = CM((M+1):(end-M),:);
		CM_e = CM((end-M+1):end,:);

		Y_b = reshape(CM_b * Y(1:w,:), [M ysn(2:end)]);

		if isempty(CM_m)
			Y_m = [];
		else
			Y_m = filter(fliplr(CM_m), 1, Y);
      Y_m = reshape(Y_m(w:end,:),[(ysn(1)-2*M) ysn(2:end)]);
    end

		Y_e = reshape(CM_e * Y((end-w+1):end,:), [M ysn(2:end)]);
    clear Y;

    Y_hat = permute([Y_b; Y_m; Y_e], order);
		clear Y_b Y_m Y_e;
	end


	if create_filter & apply_filter 
		varargout(1:2) = {Y_hat, CM};

	elseif create_filter
		varargout(1) = {CM};

	elseif apply_filter 
		varargout(1) = {Y_hat};
	end

	varargout = varargout(1:nargout);
	return;