function TestfSGolayFilt(doSpeed)
% Automatic test: fSGolayFilt
% A short test is performed, but for the main part of the test SGOLAYFILT from
% the Signal Processing Toolbox is needed.
%
% This is a routine for automatic testing. It is not needed for processing and
% can be deleted or moved to a folder, where it does not bother.
%
% TestFSGolayFilt(doSpeed)
% INPUT:
% doSpeed: If FALSE, the speed is not tested. Optional.
% OUTPUT:
% On failure the test stops with an error.
%
% Tested: Matlab 6.5, 7.7, 7.8, WinXP
% Author: Jan Simon, Heidelberg, (C) 2009-2010 matlab.THISYEAR(a)nMINUSsimon.de
% $JRev: R0p V:047 Sum:NYRQp5v70ZBk Date:08-Jun-2010 00:47:16 $
% $License: BSD (see Docs\BSD_License.txt) $
% $File: Tools\UnitTests_\TestfSGolayFilt.m $
% Initialize: ==================================================================
if nargin == 0
doSpeed = true;
end
% STARTUNITTEST creates a dialog, which allows breaking the process:
StopDlgH = getappdata(0, 'UniTest_MsgBox');
% Cannot call Matlab's SGOLAYFILT in Matlab 6 and 2010a with SINGLE array:
MatlabV = [100, 1] * sscanf(version, '%d.', 2);
SingleSupported = MatlabV > 700 && MatlabV <= 710;
% Do the work: =================================================================
disp(['==== Test fSGolayFilt ', datestr(now, 0)]);
fprintf('fSGolayCore():\n');
fSGolayCore;
fprintf('\n');
% Make noisy curves:
x(1, :) = sin(1:.1:20);
x(2, :) = cos(1:.1:20);
x = transpose(x + rand(2, length(x)) * 0.1);
y = fSGolayFilt(x, 3, 25);
% To see the differences:
%plot(x, 'r');
%hold('on');
%plot(y, 'b');
% Mean and std of the difference:
xMy = x(:) - y(:);
maxDiff = max(abs(xMy));
meanDiff = mean(xMy);
stdDiff = std(xMy);
disp('Noisy curves smoothed (desirable: max<<0.1, mean<<0.01):');
disp([' Max Diff: ', sprintf('%g', maxDiff)]);
disp([' Mean Diff: ', sprintf('%g', meanDiff)]);
disp([' Std Diff: ', sprintf('%g', stdDiff)]);
% For a comparison with SGOLAYFILT, this routine must be available:
fprintf('\n');
if isempty(which('sgolayfilt'))
disp('::: Compare fSGolayFilt with local M-version');
disp(' Speed comparisons are not really meaningful!');
sgolayfcn = @my_sgolayfilt;
cmpFuncName = 'local M-SGOLAYFILT';
else
disp('::: Compare fSGolayFilt with SGOLAYFILT of Signal-Processing-Toolbox');
sgolayfcn = @sgolayfilt;
cmpFuncName = 'Matlab''s SGOLAYFILT';
end
% Compare replies from fSGolayFilt and SGOLAYFILT ==============================
for iType = 1:2
fcn = sgolayfcn;
fcnName = cmpFuncName;
if iType == 1
disp([char(10), 'Random test data - DOUBLE:']);
x = rand(100, 10);
Tol = 1000 * eps;
else
disp([char(10), 'Random test data - SINGLE:']);
x = single(rand(100, 10));
if ~SingleSupported
fcn = @my_sgolayfilt;
fncName = 'local M-SGOLAYFILT';
end
Tol = 1.0e-6;
end
y = fSGolayFilt(x, 1, 7);
z = feval(fcn, x, 1, 7);
if isEqualTol_L(y, z, Tol)
disp([' ok: rand(100, 10), K=1, F=7 : Same result as ', fcnName]);
else
error(['TestfSGolayFilt: Results differ from ', fcnName, ':', ...
char(10), ' rand(100, 10), K=1, F=7']);
end
y = fSGolayFilt(x, 2, 5);
z = feval(fcn, x, 2, 5);
if isEqualTol_L(y, z, Tol)
disp([' ok: rand(100, 10), K=2, F=5 : Same result as ', fcnName]);
else
error(['TestfSGolayFilt: Results differ from ', fcnName, ':', ...
char(10), ' rand(100, 10), K=2, F=5']);
end
y = fSGolayFilt(x, 3, 25);
z = feval(fcn, x, 3, 25);
if isEqualTol_L(y, z, Tol)
disp([' ok: rand(100, 10), K=3, F=25: Same result as ', fcnName]);
else
error(['TestfSGolayFilt: Results differ from ', fcnName, ':', ...
char(10), ' rand(100, 10), K=3, F=25']);
end
y = fSGolayFilt(x, 5, 15);
z = feval(fcn, x, 5, 15);
if isEqualTol_L(y, z, Tol)
disp([' ok: rand(100, 10), K=5, F=15: Same result as ', fcnName]);
else
error(['TestfSGolayFilt: Results differ from ', fcnName, ':', ...
char(10), ' rand(100, 10), K=5, F=15']);
end
% Compare results from transposed vectors:
a = x(:, 1);
y1 = fSGolayFilt(a, 3, 25);
y2 = fSGolayFilt(transpose(a), 3, 25);
if isEqualTol_L(y1, transpose(y2), Tol)
disp(' ok: fSGolayFilt(transpose(X)) == transpose(fSGolayFilt(X))');
else
error('TestfSGolayFilt: Different results after transposition!');
end
% Try to test different threading strategies:
x = rand(4000, 1);
if iType == 2
x = single(x);
end
y_5_15 = fSGolayFilt(x, 5, 15);
y_7_25 = fSGolayFilt(x, 7, 25);
for Col = 1:16
xc = x(:, ones(1, Col));
yc_5_15 = fSGolayFilt(xc, 5, 15);
yc_7_25 = fSGolayFilt(xc, 7, 25);
if ~isEqualTol_L(y_5_15(:, ones(1, Col)), yc_5_15, Tol) || ...
~isEqualTol_L(y_7_25(:, ones(1, Col)), yc_7_25, Tol)
error('TestfSGolayFilt: Filtering a multi-column matrix failed');
end
end
disp(' ok: Multi-columm matrix passed the tests');
x = rand(18, 8);
if iType == 2
x = single(x);
end
if isEqualTol_L(fSGolayFilt(x, 3, 17), my_sgolayfilt(x, 3, 17), Tol)
disp(' ok: X has F+1 frames.');
else
error('TestfSGolayFilt: Failed if X has F+1 frames.');
end
x = rand(100, 4);
if iType == 2
x = single(x);
end
if isEqualTol_L(fSGolayFilt(x, 4, 5), x, Tol)
disp(' ok: X equals filtered X for K=4 and F=5.');
else
error('TestfSGolayFilt: X differs from filtered X for K=4 and F=5.');
end
if isEqualTol_L(fSGolayFilt(x, 0, 1), x, Tol)
disp(' ok: X equals filtered X for K=0 and F=1.');
else
error('TestfSGolayFilt: X differs from filtered X for K=0 and F=1.');
end
% Multi-dim X:
x = rand(11, 13, 17);
if iType == 2
x = single(x);
end
y = fSGolayFilt(x, 3, 7, [], 1);
y_ = fSGolayFilt(x, 3, 7, [], []);
y2 = my_sgolayfilt(x, 3, 7);
if ~isEqualTol_L(y, y2, Tol) || ~isequal(y_, y)
error('TestfSGolayFilt: Dim=1 on multi-dim X failed');
end
y = fSGolayFilt(x, 3, 7, [], 2);
yt = my_sgolayfilt(permute(x, [2, 1, 3]), 3, 7);
y2 = permute(reshape(yt, [13, 11, 17]), [2, 1, 3]);
if ~isEqualTol_L(y, y2, Tol)
error('TestfSGolayFilt: Dim=2 on multi-dim X failed');
end
y = fSGolayFilt(x, 3, 7, [], 3);
yt = my_sgolayfilt(permute(x, [3, 2, 1]), 3, 7);
y2 = permute(reshape(yt, [17, 13, 11]), [3, 2, 1]);
if ~isEqualTol_L(y, y2, Tol)
error('TestfSGolayFilt: Dim=3 on multi-dim X failed');
end
disp(' ok: All 3 dims of a 3D array filtered correctly.');
end
% Speed: -----------------------------------------------------------------------
disp([char(10), '== Test speed:']);
minT = 10 * eps; % Avoid division by zero
if doSpeed
% Number of loops adjusted to processor speed:
x = rand(100, 10);
iLoop = 0;
startTime = cputime;
while cputime - startTime < 1.0
z = fSGolayFilt(x, 5, 15); %#ok<*NASGU>
clear('z'); % Reduce influence of re-using memory
iLoop = iLoop + 1;
end
nLoops = 100 * ceil(iLoop / ((cputime - startTime) * 200));
else
% At least 2 times to give the MEX files a chance to crash...
disp(' Use at least 2 loops (Displayed times are random!)');
nLoops = 2;
end
K = 3;
F = 25;
fprintf([' Fixed filter parameters: K=%d, F=%d, time in [sec], ', ...
'%d loops\n'], K, F, nLoops);
avSpeed = [];
for C = [1, 2, 4, 8, 16] % Number of signals
for L = [100, 200, 400, 1600] % Length of signals
x = rand(L, C);
if L == 1600
factor = 16;
else
factor = 4;
end
tic;
for i = 1:nLoops / factor
z = feval(sgolayfcn, x, K, F); clear('z');
end
sgt = minT + toc * factor;
tic;
for i = 1:nLoops
fz = fSGolayFilt(x, K, F); clear('fz');
end
ft = minT + toc;
% Compare the results again, who knows:
if ~isEqualTol_L(sgolayfilt(x, K, F), fSGolayFilt(x, K, F), 1000*eps)
disp('*** WARNING: Different results!');
end
disp([' X=[', sprintf('%4d %2d', L, C), ']: ', ...
cmpFuncName, ': ', sprintf('%5.2f', sgt), ...
' fSGolayFilt: ', sprintf('%5.2f', ft), ...
' ==> ', sprintf('%4.1f', 100 * ft / sgt), '%']);
avSpeed(end + 1) = 100 * ft / sgt; %#ok<AGROW>
% Break process on demand:
drawnow;
if ~all(ishandle(StopDlgH))
disp([mfilename, ': User aborted process.']);
return;
end
end
end
fprintf(' Average: %.2f +- %.2f %% of %s\n', ...
mean(avSpeed), std(avSpeed), cmpFuncName);
disp([char(10), ' Fixed X=[200, 8]']);
x = rand(200, 8);
avSpeed = [];
for K = [3, 5, 7]
for F = [9, 15, 25]
tic;
for i = 1:nLoops / 2
z = feval(sgolayfcn, x, K, F); clear('z');
end
sgt = minT + toc * 2;
tic;
for i = 1:nLoops
fz = fSGolayFilt(x, K, F); clear('fz');
end
ft = minT + toc;
% Compare the results again, who knows:
if ~isEqualTol_L(sgolayfilt(x, K, F), fSGolayFilt(x, K, F), 1000*eps)
disp('*** WARNING: Different results!');
end
disp([' K=', sprintf('%d', K), ', F=', sprintf('%2d', F), ...
': ', cmpFuncName, ': ', sprintf('%5.2f', sgt), ...
' fSGolayFilt: ', sprintf('%5.2f', ft), ...
' ==> ', sprintf('%4.1f', 100 * ft / sgt) '%']);
avSpeed(end+1) = 100 * ft / sgt; %#ok<AGROW>
% Break process on demand:
drawnow;
if ~all(ishandle(StopDlgH))
disp([mfilename, ': User aborted process.']);
return;
end
end
end
fprintf(' Average: %.2f +- %.2f %% of %s\n\n', ...
mean(avSpeed), std(avSpeed), cmpFuncName);
% Number of loops adjusted to processor speed: ---------------------------------
for iType = 1:2
for k = [1, 2, 4]
if SingleSupported
x = rand(1E6, k);
sizeStr = '1E6';
else
x = rand(1E5, k);
sizeStr = '1E5';
end
if doSpeed
iLoop = 0;
startTime = cputime;
while cputime - startTime < 1.0
z = fSGolayFilt(x, 3, 25);
clear('z'); % Reduce influence of JIT
iLoop = iLoop + 1;
end
nLoops = ceil(iLoop / (cputime - startTime));
else
nLoops = 2;
end
timeFactor = 1;
if iType == 1
x = double(x);
fprintf(' Large X=double(rand(%s, %d)), %d loops\n', ...
sizeStr, k, nLoops);
fcn = sgolayfcn;
fcnName = cmpFuncName;
else
fprintf(' Large X=single(rand(%s, %d)), %d loops\n', ...
sizeStr, k, nLoops);
x = single(x);
if ~SingleSupported
fprintf(' (SINGLEs are no supported in this Matlab version\n');
fcn = @my_sgolayfilt;
fcnName = 'local M-SGOLAYFILT';
timeFactor = 10;
end
end
for iL = 1:3
switch iL
case 1
K = 3;
F = 25;
case 2
K = 5;
F = 15;
case 3
K = 4;
F = 9;
otherwise
error('Programming error: Bad SWITCH');
end
tic;
for i = 1:nLoops
fz = fSGolayFilt(x, K, F); clear('fz');
end
ft = minT + toc;
tic;
for i = 1:nLoops / timeFactor
z = feval(fcn, x, K, F); clear('z');
end
sgt = minT + toc * timeFactor;
% Compare the results:
y1 = feval(fcn, x, K, F);
y2 = fSGolayFilt(x, K, F);
if ~isEqualTol_L(y1, y2, 100 * sqrt(eps))
disp('*** WARNING: Different results!');
end
disp([' K=', sprintf('%d', K), ', F=', sprintf('%2d', F), ...
': ', fcnName, ': ', sprintf('%5.2f', sgt), ...
' fSGolayFilt: ', sprintf('%5.2f', ft), ...
' ==> ', sprintf('%4.1f', 100 * ft / sgt) '%']);
end
end
end
% Goodbye:
fprintf('\nfSGolayFilt seems to work well.\n');
return;
% ******************************************************************************
function Equal = isEqualTol_L(x, y, Tol)
% A comparison with some tolerance
if max(abs(double(x(:)) - double(y(:)))) <= Tol
Equal = true;
else
Equal = false;
end
return;
% ******************************************************************************
function Y = my_sgolayfilt(X, K, F, W)
% Savitzky-Golay polynomial filtering [No MEX]
% A local version, faster than SGOLAYFILT of Matlab 5.3
% X can be DOUBLE or SINGLE.
% MatLab 5.3, 6.0, matlab.THISYEAR(a)nMINUSsimon.de, (C) 2003
% Was: JRev: R0a Build004 Sum:137B3D Date:26-Feb-2003 00:09
% Initialize: ==================================================================
% Global Interface: ------------------------------------------------------------
FuncName = 'my_sgolayfilt';
% Initial values: --------------------------------------------------------------
[sX1, sX2] = size(X);
if sX1 == 1 % Make X a column oriented matrix if it is a vector:
X = X(:);
sX1 = sX2;
sX2 = 1;
end
% Program Interface: -----------------------------------------------------------
error(nargchk(3, 4, nargin));
% Test input arguments:
if round(F) ~= F || rem(F, 2) ~= 1
error([FuncName ': Frame length F must be an odd integer.']);
end
if round(K) ~= K
error([FuncName ': Polynomial order K must be an integer.']);
end
if K > F-1
error([FuncName ': Order K must be smaller than frame length F.']);
end
if F > sX1
error([FuncName ': Input signal must be longer than frame length F.']);
end
% No SINGLE arithmetics in Matlab 6.5:
XtoDouble = isa(X, 'single') && sscanf(version, '%f', 1) < 7.0;
if XtoDouble
X = double(X);
Y(sX1, sX2) = single(0);
else
Y = zeros(sX1, sX2); % Allocate memory for output matrix
end
% Do the work: =================================================================
if nargin < 4
% No weighting matrix:
% Compute the projection matrix B (faster SGOLAY without weight)
S = localVander(F, K + 1); % Vandermonde like matrix
R = chol(transpose(S) * S); % Cholesky factor
A = S / R;
B = A * transpose(A);
else % nargin == 4:
if length(W) ~= F
error([FuncName ': Weight vector and F must have same length.']);
end
if min(W) <= 0
error([FuncName ': Weight vector must have positive elements.']);
end
% Projection matrix B:
W = diag(W);
S = localVander(F, K + 1); % Vandermonde like matrix
R = chol(transpose(S) * W * S); % Cholesky factor
A = S / R;
B = A * transpose(A) * W;
end
% Compute smoothing result for each column: ------------------------------------
sF = sX1 - F - 1;
dF = (F - 1) / 2 + 1;
% Initial and final smoothed data:
Y(1:dF, :) = B(end:-1:dF, :) * X(F:-1:1, :);
Y(sX1 - dF + 1:sX1, :) = B(dF:-1:1, :) * X(sX1:-1:sX1 - F + 1, :);
% xBUF contains the moving windows with F frames:
BUF = [2:sX1 - F; ones(F-1, sF)]; % Matrix of indices
xBUF = X(cumsum(BUF, 1), :); % Big 2D matrix of windowed values
% Compute the steady state output:
Y(dF + 1:sX1 - dF, :) = reshape(B(dF, :) * reshape(xBUF, F, sF * sX2), sF, sX2);
return;
% ******************************************************************************
function V = localVander(F, N)
% Special Vandermonde matrix (see VANDER)
L = (F-1) / 2;
c = transpose(-L:L);
V = ones(F, N);
P = ones(F, 1);
for i = 2:N
P = P .* c;
V(:, i) = P;
end
return;
