fSGolayFilt(X, K, F, W, Dim)

function Y = fSGolayFilt(X, K, F, W, Dim)
% Fast Savitzky-Golay polynomial filtering [Call MEX]
% Y = fSGolayFilt(X, K, F, [W], [Dim])
% The uniformly spaced signal X is smoothed by computing a polynomial of order K
% over windows with F frames in least-squares distance to the original signal.
% This method was developped by Savitzky and Golay.
% INPUT:
%   X: Vector [1 x nFrame], [nFrame x 1] or matrix [nFrame x nChannel].
%      Class: SINGLE or DOUBLE.
%   K: Order of the filter polynomial as integer.
%   F: Odd number of frames of filter window. F > K and F <= nFrame.
%   W: Real positive weight factors for the least-squares method.
%      Optional, default: [].
%   Dim: Dimension to operate on. Optional, default: first non-singelton
%      dimension.
%
% OUTPUT:
%   Y: Array of filtered values with same size as input X. No phase distortion
%      appears. Smoothing is stronger for smaller K and larger F. For K = F-1, 
%      the output equals the input.
%
% Input and output of this function equal the function SGOLAYFILT of the
% Signal-Processing-Toolbox written by R. Losada. Under Matlab 2009a this C-Mex
% implementation is about 85% faster for short [100 x 1] vectors and 30% faster
% for [1E7 x 1] vectors.
% The last used projection matrix stored persistently to accelerate the repeated
% application of a specific set of filter parameters to different data.
%
% EXAMPLES:
%   X = sin(0:0.1:10 * pi);       % sin wave
%   Y = X + 0.1 * rand(size(X));  % add noise
%   Z = fSGolayFilt(XX, 3, 25);   % smooth, polynomials of deg 3 over 25 points
%   plot(1:length(Y), Y, 'b', 1:length(Z), Z, 'r');
%
% MULTI-THREADING:
%   At compile time the largest number of cores is defined by the directive
%   MAX_NTHREAD (default: 8). For the current computer, the number of threads is
%   limited by the environment variable OMP_NUM_THREADS (fallback for Matlab
%   6.5: NUMBER_OF_PROCESSORS). The user can limit the number of threads
%   manually 
%   The compiled Mex function starts mutliple threads to solve the problem. The 
%   optimal number of threads is estimated automatically using the size of the
%   input, the window length F and the parameter THREAD_CONST. The largest
%   number of threads MAX_NTHREAD must be defined at compile time. 
%
% COMPILE:
%   See "fSGolayCore.c" and "fSGolayCore_ST.c".
%   The LCC or OWC compilers do not create suffciently fast programs for this
%   C-Mex.
%
% TEST: Run TestfSGolayFilt after compiling to check validity and speed.
%   This test function is not needed for running fSGolayFilt.
%
% See also SGOLAYFILT, SGOLAY, VANDER, MEDFILT1, FILTER.
%
% Tested: Matlab 6.5, 7.7, 7.8, WinXP
%         64 bit compatibility is assumed.
%         Single threaded: Windows, Linux, OSX. Multi threaded: Windows.
%         Compilers:  BCC5.5, LCC2.4/3.8, Open Watcom 1.8, MSVC++ 2008
% Author: Jan Simon, Heidelberg, (C) 2006-2010 matlab.THISYEAR(a)nMINUSsimon.de

% REFERENCES:
% [1] Sophocles J. Orfanidis: "Introduction to Signal Processing"
%     Prentice-Hall, 1995, Chapter 8.
% [2] 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.
% [3] K. L. Ratzlaff and J. T. Johnson: "Computation of Two-Dimensional
%     Polynomial Least-Squares Convolution Smoothing Integers",
%     Analytical Chemistry, vol. 61, pp. 1303-1305, 1989.
% [4] J. E. Kuo, H. Wang, and S. Pickup,
%     "Multidimensional Least-Squares Smoothing Using Orthogonal Polynomials"
%     Analytical Chemistry, vol. 63, pp. 630-635, 1991.

% $JRev: R0s V:063 Sum:q6hRM5cB56s4 Date:08-Jun-2010 00:47:16 $
% $License: BSD $
% $UnitTest: TestfSGolayFilt $
% $File: Tools\GLMath\fSGolayFilt.m $
% History:
% 043: 28-Oct-2006 00:07, Faster Vander_L.
% 047: 27-May-2008 10:43, Faster CAT of output, remember last [K, F, B].
%      The remembering of the last matrix B is not limited to the parameters K=3
%      and F=25 anymore. Handle empty input.
% 054: 04-Mar-2010 23:42, Initial and final phase in Mex also => 10% faster.
%      The overhead of the main function in the M-file is just about 1% for real
%      world problems.
% 062: 07-Jun-2010 01:10, Multi-threading on Windows, input [Dim], N-D input.
%      B is DOUBLE even when X is a single. S*inv(R) => S/R.

% Initialize: ==================================================================
% Global Interface: ------------------------------------------------------------
% Container for projection matrix B for last used parameters:
persistent Last_KF Last_B

% Initial values: --------------------------------------------------------------
if isempty(X)
   Y = X;
   return;
end

% Program Interface: -----------------------------------------------------------
nArg = nargin;  % Save time for repeated calls of NARGIN
if or(nArg < 3, nArg > 5)
    error(['*** ', mfilename, ': Bad number of inputs.']);
end

% Check inputs:
if rem(F, 2) ~= 1  % | round(F) ~= F  ( rem(X,2)==1 only for integers! )
   error(['*** ', mfilename, ': Window length F must be an odd integer.']);
end
if round(K) ~= K
   error(['*** ', mfilename, ': Polynomial order K must be an integer.']);
end
if K > F - 1
   error(['*** ', mfilename, ...
         ': Order K must be smaller than window length F.']);
end

% Default values for optional arguments W and Dim:
if nArg < 4
   W   = [];
   Dim = [];
elseif nArg < 5
   Dim = [];
end

% Put Dim to the first dimension:
sizeX = size(X);
if ndims(X) == 2    % Standard case: vector or matrix
   if isempty(Dim)
      doPermute = false;
      doReshape = (sizeX(1) == 1);
      if doReshape
         X = X(:);
      end
   elseif Dim > 2
      error(['*** ', mfilename, ': [Dim] out of range.']);
   else
      doReshape = false;
      doPermute = (Dim == 2);
      if doPermute
         perm = [2, 1];
         X    = transpose(X);
      end
   end
else                % More than 2 dimensions:
   doPermute = false;
   doReshape = true;
   if isempty(Dim)  % First non-singleton dimension:
      non1Ind = find(sizeX ~= 1);
      if length(non1Ind) && non1Ind(1) > 1
         tmp = [sizeX, 1];
         X   = reshape(X, tmp(non1Ind(1):length(tmp)));
      end
   elseif Dim > 1   % Dim specified:
      if Dim > ndims(X)
         error(['*** ', mfilename, ': [Dim] out of range.']);
      end
      doPermute = true;
      perm      = [Dim, 1:Dim - 1, Dim + 1:ndims(X)];
      X         = permute(X, perm);
      sizeX     = sizeX(perm);
   end
end

if F >= size(X, 1)
   error(['*** ', mfilename, ': Signal must be longer than window length F.']);
end

% Do the work: =================================================================
if isempty(W)                      % No weighting matrix:
   KF = [K, F];
   if isequal(KF, Last_KF)
      B = Last_B;                  % Same parameter as last call
   else                            % Calculate matrix B:
      S = Vander_L(F, K + 1);      % Vandermonde like matrix
      R = chol(transpose(S) * S);  % Cholesky factor
      A = S / R;
      B = A * transpose(A);
      
      Last_B  = B;                 % Remember this B for the future
      Last_KF = KF;
   end
else                               % Weight matrix defined:
   if length(W) ~= F
      error(['*** ', mfilename, ': Weight vector needs length of F.']);
   end
   if any(W <= 0)
      error(['*** ', mfilename, ': Weight vector needs positive elements.']);
   end
   
   % Projection matrix B:
   W = diag(W);
   S = Vander_L(F, K + 1);         % Vandermonde like matrix
   R = chol(transpose(S) * W * S); % Cholesky factor
   A = S / R;
   B = A * transpose(A) * W;
end

% Let the MEX perform the calculations:
Y = fSGolayCore(X, F, B);
% [Y, n] = fSGolayCore(X, F, B);  % Replies the number of used threads

% Reconstruct original vector dimension:
if doReshape
   Y = reshape(Y, sizeX);
end
if doPermute
   Y = ipermute(Y, perm);
end

return;

% ******************************************************************************
function V = Vander_L(F, N)
% Special Vandermonde matrix (see VANDER)
% Just the elements needed for the Savitzky-Golay-Filter are computed:
L = (F - 1) / 2;
c = reshape(-L:L, [], 1);
P = c(:, ones(1, N));
P(:, 1) = 1.0;
V = cumprod(P, 2);

return;