Code covered by the BSD License  

Highlights from
ButterParam

ButterParam

by

 

31 Jan 2011 (Updated )

Store persistent list of BUTTER parameters for faster access

uTest_ButterParam(doSpeed)
function uTest_ButterParam(doSpeed)
% Unit test: ButterParam
% 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.
%
% uTest_ButterParam(doSpeed)
% INPUT:
%   doSpeed: If ANY(doSpeed) is TRUE, the speed is measured.
% OUTPUT:
%   On failure the test stops with an error.
%
% Tested: Matlab 6.5, 7.7, 7.8, WinXP
% Author: Jan Simon, Heidelberg, (C) 2009-2011 matlab.THISYEAR(a)nMINUSsimon.de

% $JRev: R-k V:010 Sum:NrHnk2zXMIcW Date:12-Feb-2011 00:46:28 $
% $License: NOT_RELEASED $
% $File: Tools\UnitTests_\uTest_ButterParam.m $

% Initialize: ==================================================================
% Global Interface: ------------------------------------------------------------
ErrID = ['JSimon:', mfilename];

% Initial values: --------------------------------------------------------------
% Program Interface: -----------------------------------------------------------
if nargin == 0
   doSpeed = true;
end

if doSpeed
   TestTime = 1.0;  % Time for random tests
else
   TestTime = 0.1;
end

% User Interface: --------------------------------------------------------------
% Do the work: =================================================================
disp(['==== Test ButterParam  ', datestr(now, 0)]);
disp(['  Version: ', which('ButterParam')]);
pause(0.01);  % Flush events

% It is not absolutely necessary to have the Signal-Processing-Toolbox to use
% this function, as long as all used parameters are store in the preferences
% file of ButterParam. But the SPT command BUTTER is needed to create the
% preferences at first.
if isempty(which('butter'))
   disp('! Signal-Processing-Toolbox not available !');
   disp('! Testing is not possible !');
   return;
end

% Consider slightly different values in Matlab 6.5 and 2009a:
Tol = 100 * eps;

% Known answer tests: ----------------------------------------------------------
% Order 1:
[B, A] = ButterParam(1, 0.5);
if isEqualTol(B, [0.5, 0.5], Tol) && isEqualTol(A, [1.0, 0.0], Tol)
   disp('  ok: (1, 0.5)');
else
   error([ErrID, 'KATfail'], 'Failed for (1, 0.5)');
end

[B, A] = ButterParam(1, 0.5, 'low');
if isEqualTol(B, [0.5, 0.5], Tol) && isEqualTol(A, [1.0, 0.0], Tol)
   disp('  ok: (1, 0.5, low)');
else
   error([ErrID, 'KATfail'], 'Failed for (1, 0.5, low)');
end

[B, A] = ButterParam(1, 0.5, 'high');
if isEqualTol(B, [0.5, -0.5], Tol) && isEqualTol(A, [1.0, 0.0], Tol)
   disp('  ok: (1, 0.5, high)');
else
   error([ErrID, 'KATfail'], 'Failed for (1, 0.5, high)');
end

[B, A] = ButterParam(1, 0.3);
Bw     = [0.337540151883547, 0.337540151883547];
Aw     = [1, -0.324919696232906];
if isEqualTol(B, Bw, Tol) && isEqualTol(A, Aw, Tol)
   disp('  ok: (1, 0.3)');
else
   error([ErrID, 'KATfail'], 'Failed for (1, 0.3)');
end

[B, A] = ButterParam(1, 0.3, 'low');
if isEqualTol(B, Bw, Tol) && isEqualTol(A, Aw, Tol)
   disp('  ok: (1, 0.3, low)');
else
   error([ErrID, 'KATfail'], 'Failed for (1, 0.3, low)');
end

[B, A] = ButterParam(1, 0.3, 'high');
Bw     = [0.6624598481164532, -0.6624598481164532];
Aw     = [1, -0.3249196962329063];
if isEqualTol(B, Bw, Tol) && isEqualTol(A, Aw, Tol)
   disp('  ok: (1, 0.3, high)');
else
   error([ErrID, 'KATfail'], 'Failed for (1, 0.3, high)');
end

[B, A] = ButterParam(1, [0.4, 0.6], 'bandpass');
Bw     = [0.2452372752527857, 0, -0.2452372752527857];
Aw     = [1, 0, 0.5095254494944287];
if isEqualTol(B, Bw, Tol) && isEqualTol(A, Aw, Tol)
   disp('  ok: (1, [0.4, 0.6], band)');
else
   error([ErrID, 'KATfail'], 'Failed for (1, [0.4, 0.6], band)');
end

[B, A] = ButterParam(1, [0.4, 0.6], 'stop');
Bw     = [0.7547627247472142, 0, 0.7547627247472142];
Aw     = [1, 0, 0.5095254494944287];
if isEqualTol(B, Bw, Tol) && isEqualTol(A, Aw, Tol)
   disp('  ok: (1, [0.4, 0.6], stop)');
else
   error([ErrID, 'KATfail'], 'Failed for (1, [0.4, 0.6], stop)');
end

% 2nd order: -------------------------------------------------------------------
[B, A] = ButterParam(2, 0.5);
Bw     = [0.2928932188134524, 0.5857864376269049, 0.2928932188134524];
Aw     = [1, 0, 0.1715728752538099];
if isEqualTol(B, Bw, Tol) && isEqualTol(A, Aw, Tol)
   disp('  ok: (2, 0.5)');
else
   error([ErrID, 'KATfail'], 'Failed for (2, 0.5)');
end

[B, A] = ButterParam(2, 0.5, 'low');
if isEqualTol(B, Bw, Tol) && isEqualTol(A, Aw, Tol)
   disp('  ok: (2, 0.5, low)');
else
   error([ErrID, 'KATfail'], 'Failed for (2, 0.5, low)');
end

[B, A] = ButterParam(2, 0.5, 'high');
Bw     = [0.2928932188134526, -0.5857864376269052, 0.2928932188134526];
Aw     = [1, 0, 0.1715728752538099];
if isEqualTol(B, Bw, Tol) && isEqualTol(A, Aw, Tol)
   disp('  ok: (2, 0.5, high)');
else
   error([ErrID, 'KATfail'], 'Failed for (2, 0.5, high)');
end

[B, A] = ButterParam(2, 0.3);
Bw     = [0.131106439916626, 0.262212879833252, 0.131106439916626];
Aw     = [1, -0.747789178258503, 0.272214937925007];
if isEqualTol(B, Bw, Tol) && isEqualTol(A, Aw, Tol)
   disp('  ok: (2, 0.3)');
else
   error([ErrID, 'KATfail'], 'Failed for (2, 0.3)');
end

[B, A] = ButterParam(2, 0.3, 'low');
if isEqualTol(B, Bw, Tol) && isEqualTol(A, Aw, Tol)
   disp('  ok: (2, 0.3, low)');
else
   error([ErrID, 'KATfail'], 'Failed for (2, 0.3, low)');
end

[B, A] = ButterParam(2, 0.3, 'high');
Bw     = [0.5050010290458777, -1.010002058091756, 0.5050010290458777];
Aw     = [1, -0.7477891782585036, 0.2722149379250074];
if isEqualTol(B, Bw, Tol) && isEqualTol(A, Aw, Tol)
   disp('  ok: (2, 0.3, high)');
else
   error([ErrID, 'KATfail'], 'Failed for (2, 0.3, high)');
end

[B, A] = ButterParam(2, [0.4, 0.6], 'bandpass');
Bw     = [0.06745527388907176, 0, -0.1349105477781435, 0, 0.06745527388907176];
Aw     = [1, 0, 1.142980502539903, 0, 0.4128015980961897];
% Aw     = [1, 0, 1.1429805025399, 0, 0.412801598096188];  % Matlab 2009a
if isEqualTol(B, Bw, Tol) && isEqualTol(A, Aw, Tol)
   disp('  ok: (2, [0.4, 0.6], band)');
else
   error([ErrID, 'KATfail'], 'Failed for (2, [0.4, 0.6], band)');
end

[B, A] = ButterParam(2, [0.4, 0.6], 'stop');
Bw     = [0.6389455251590226, 0, 1.277891050318045, 0, 0.6389455251590226];
Aw     = [1, 0, 1.142980502539901, 0, 0.4128015980961889];
if isEqualTol(B, Bw, Tol) && isEqualTol(A, Aw, Tol)
   disp('  ok: (2, [0.4, 0.6], stop)');
else
   error([ErrID, 'KATfail'], 'Failed for (2, [0.4, 0.6], stop)');
end

% 3rd order: -------------------------------------------------------------------
[B, A] = ButterParam(3, 0.5);
Bw     = [0.1666666666666666, 0.4999999999999999, 0.4999999999999999, ...
   0.1666666666666666];
Aw     = [1, 0, 0.3333333333333334, 0];
if isEqualTol(B, Bw, Tol) && isEqualTol(A, Aw, Tol)
   disp('  ok: (3, 0.5)');
else
   error([ErrID, 'KATfail'], 'Failed for (3, 0.5)');
end

[B, A] = ButterParam(3, 0.5, 'low');
if isEqualTol(B, Bw, Tol) && isEqualTol(A, Aw, Tol)
   disp('  ok: (3, 0.5, low)');
else
   error([ErrID, 'KATfail'], 'Failed for (3, 0.5, low)');
end

[B, A] = ButterParam(3, 0.5, 'high');
Bw     = [0.1666666666666667, -0.5000000000000001, 0.5000000000000001, ...
   -0.1666666666666667];
Aw     = [1, 0, 0.3333333333333334, 0];
if isEqualTol(B, Bw, Tol) && isEqualTol(A, Aw, Tol)
   disp('  ok: (3, 0.5, high)');
else
   error([ErrID, 'KATfail'], 'Failed for (3, 0.5, high)');
end

[B, A] = ButterParam(3, 0.3);
Bw     = [0.0495329963572532, 0.1485989890717596, 0.1485989890717596, ...
   0.0495329963572532];
Aw     = [1, -1.161917483671733, 0.6959427557896512, -0.1377613012598929];
if isEqualTol(B, Bw, Tol) && isEqualTol(A, Aw, Tol)
   disp('  ok: (3, 0.3)');
else
   error([ErrID, 'KATfail'], 'Failed for (3, 0.3)');
end

[B, A] = ButterParam(3, 0.3, 'low');
if isEqualTol(B, Bw, Tol) && isEqualTol(A, Aw, Tol)
   disp('  ok: (3, 0.3, low)');
else
   error([ErrID, 'KATfail'], 'Failed for (3, 0.3, low)');
end

[B, A] = ButterParam(3, 0.3, 'high');
Bw     = [0.3744526925901596, -1.123358077770479, 1.123358077770479, ...
   -0.3744526925901596];
Aw     = [1, -1.161917483671733, 0.695942755789651, -0.1377613012598929];
if isEqualTol(B, Bw, Tol) && isEqualTol(A, Aw, Tol)
   disp('  ok: (3, 0.3, high)');
else
   error([ErrID, 'KATfail'], 'Failed for (3, 0.3, high)');
end

[B, A] = ButterParam(3, [0.4, 0.6], 'bandpass');
Bw     = [0.01809893300751436, 0, -0.05429679902254309, 0, ...
      0.05429679902254309, 0, -0.01809893300751436];
Aw     = [1, 0, 1.760041880343171, 0, 1.182893262037834, 0, 0.2780599176345473];
% Aw     = [1, 0, 1.1429805025399, 0, 0.412801598096188];  % Matlab 2009a
if isEqualTol(B, Bw, Tol) && isEqualTol(A, Aw, Tol)
   disp('  ok: (3, [0.4, 0.6], band)');
else
   error([ErrID, 'KATfail'], 'Failed for (3, [0.4, 0.6], band)');
end

[B, A] = ButterParam(3, [0.4, 0.6], 'stop');
Bw     = [0.5276243825019429, 0, 1.582873147505829, 0, 1.582873147505829, ...
   0, 0.5276243825019429];
Aw     = [1, 0, 1.760041880343167, 0, 1.182893262037828, 0, 0.2780599176345457];
if isEqualTol(B, Bw, Tol) && isEqualTol(A, Aw, Tol)
   disp('  ok: (3, [0.4, 0.6], stop)');
else
   error([ErrID, 'KATfail'], 'Failed for (3, [0.4, 0.6], stop)');
end

% Now check an entry appearing not at the end of the list again:
[B, A] = ButterParam(2, 0.3);
Bw     = [0.131106439916626, 0.262212879833252, 0.131106439916626];
Aw     = [1, -0.747789178258503, 0.272214937925007];
if isEqualTol(B, Bw, Tol) && isEqualTol(A, Aw, Tol)
   disp('  ok: (2, 0.3) repeated');
else
   error([ErrID, 'KATfail'], 'Failed for (2, 0.3) repeated');
end

% Bad input test: --------------------------------------------------------------
fprintf('\n== Bad input tests (provoked errors):\n');
tooLazy = false;
try
   [B, A]  = ButterParam(1);
   tooLazy = true;
catch
   fprintf('  ok: 1 input caught: %s\n', lasterr);
end
if tooLazy
   error([ErrID, 'TooLazy'], 'Bad data accepted: ButterParam(1)');
end

try
   [B, A]  = ButterParam('wrongCommand');
   tooLazy = true;
catch
   fprintf('  ok: wrongCommand caught: %s\n', lasterr);
end
if tooLazy
   error([ErrID, 'TooLazy'], 'Unknown command accepted');
end

try
   [B, A]  = ButterParam(1, 1.1);
   tooLazy = true;
catch
   fprintf('  ok: Bad frequency caught: %s\n', lasterr);
end
if tooLazy
   error([ErrID, 'TooLazy'], 'Bad frequency accepted');
end

% No analog filter creation:
try
   [B, A]  = ButterParam(1, 0.5, 'low', 's');
   tooLazy = true;
catch
   fprintf('  ok: Bad frequency caught: %s\n', lasterr);
end
if tooLazy
   error([ErrID, 'TooLazy'], '4 inputs accepted');
end

% Speed test: ------------------------------------------------------------------
fprintf('\n== Speed tests:\n');

disp('Relation of times for getting paramters and the filtering:');
for DataLen = [100, 1000, 10000, 100000]
   Data = rand(DataLen, 1);
   
   butterTime = 0;
   paramTime  = 0;
   filterTime = 0;
   
   endTime = cputime + 2.0;
   while cputime < endTime
      tic;
      [B, A]    = ButterParam(3, 0.6, 'low');
      paramTime = paramTime + toc;
      clear('A', 'B');
      
      tic;
      [B, A]     = butter(3, 0.6, 'low');  %#ok<*NASGU>
      butterTime = butterTime + toc;
      
      tic;
      Data2      = filter(B, A, Data);
      filterTime = filterTime + toc;
      
      clear('A', 'B', 'Data2');
   end
   
   allButterTime = filterTime + butterTime;
   fprintf('  Data=[%6d x 1]: %2d%% BUTTER,      %2d%% FILTER\n', ...
      DataLen, ...
      round(100 * butterTime / allButterTime), ...
      round(100 * filterTime / allButterTime));

   allParamTime = filterTime + paramTime;
   fprintf('                     %2d%% ButterParam, %2d%% FILTER\n', ...
      round(100 * paramTime  / allParamTime), ...
      round(100 * filterTime / allParamTime));
end

fprintf('\nGet parameters:\n');
Pass = 'low';
W    = 0.5;
for iOrder = [1, 2, 4, 8, 16]
   fprintf('N=%d, W=%g, P=%s\n', iOrder, W, Pass);
      
   loops = 0;
   tic;
   elapsed = toc;
   while elapsed < TestTime
      [B, A]  = butter(iOrder, W, Pass);  %#ok<*NASGU>
      clear('A', 'B');
      loops   = loops + 1;
      elapsed = toc;
   end
   fprintf('  BUTTER:      %5d loops/s\n', round(loops / TestTime));
   
   loops = 0;
   tic;
   elapsed = toc;
   while elapsed < TestTime
      [B, A]  = ButterParam(iOrder, W, Pass);
      clear('A', 'B');
      loops   = loops + 1;
      elapsed = toc;
   end
   fprintf('  ButterParam: %5d loops/s\n', round(loops / TestTime));
end

% Ready:
fprintf('\n== ButterParam passed the test\n');

return;

% ******************************************************************************
% Dummy function to shadow ButterParam by the original BUTTEr function:
% function [B,A] = ButterParam(varargin)
% [B, A] = butter(varargin{:});

% ******************************************************************************
function Equal = isEqualTol(x, y, Tol)
% Tested: Matlab 6.5, 7.7, 7.8, WinXP
% Author: Jan Simon, Heidelberg, (C) 2007-2011 matlab.THISYEAR(a)nMINUSsimon.de

Equal = false;
if isequal(size(x), size(y))
   xMy = double(x) - double(y);   % No operations on SINGLEs in Matlab 6!
   
   % Same as "if all(abs(xMy(:)) <= Tol)", but faster:
   if all(or((abs(xMy) <= Tol), (x == y)))   % is FALSE for NaNs
       Equal = true;
   end
end

return;

Contact us