source: hypertools/hypertools/sg.m

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;