CIC decimation filter design model for Matlab/Octave (computing register pruning)

cic_decimator_model.m

% Matlab script modeling CIC decimation filter design (improved)
% as described in article
%
% E. B. Hogenauer
% "An economical class of digital filters for decimation and interpolation"
% IEEE Transactions on Acoustics, Speech and Signal Processing,
% ASSP29(2):155-162, 1981.
%
% Improvements introduce correction of negative Bj values and calculating
% extra properties.
%
% Script compatible with both Matlab and GNU Octave.
% It also shows filter analysis when available
% (for Matlab - 'Signal Processing Toolbox' Add-On installed,
%  for GNU Octave - 'signal' package installed and loaded).
%
% For Matlab users.
% Looks like script core functionality is already implemented in
% dsp.CICDecimator system object (available with 'DSP System Toolbox'
% and 'Fixed-Point Designer' Add-Ons installed). At least, I checked that
% it determines same stage parameters as this script with all default
% parameters except overrideB = [].
% To reproduce this behavior:
% - create object "cicDecim = dsp.CICDecimator(R, M, N)";
% - set "cicDecim.FixedPointDataType = 'Minimum section word lengths'";
% - set "cicDecim.OutputWordLength = Bout";
% - construct fixed-point type "nt = numerictype(true, Bin, 0)";
% - call "[WLs, FLs] = getFixedPointInfo(cicDecim, nt)";
% - WLs should be equal to accum widths;
% - FLs should be equal to -B.
% The only value this scripts adds to it is calculating error parameters
% (μᴛ, σᴛ, error RMS), but these calculations may be re-implemented using
% higher-level 'numerictype', 'quantizer', 'errmean' and 'errvar' functions.
%
% For GNU Octave users.
% If script fails with error on undefined 'table' and 'prettyprint'
% functions, then install and load 'tablicious' package.

isoctave = (exist('OCTAVE_VERSION', 'builtin') ~= 0);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input parameters definition (with default values)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

assertintegerscalar = @(v) assert(   isscalar(v) ...
                                  && isnumeric(v) ...
                                  && (v == round(v)), ...
                                  [inputname(1) ' is valid integer']);

% Rate change factor
if exist('R', 'var')
  assertintegerscalar(R);
  assert((R > 1), "R > 1");
else
  % default value from 'Design example' in article
  R = 25;
end

% Differential delay
if exist('M', 'var')
  assertintegerscalar(M);
  assert((M > 0), "M > 0");
else
  % default value from 'Design example' in article
  M = 1;
end

% Number of stages in each filter section (order)
if exist('N', 'var')
  assertintegerscalar(N);
  assert((N > 0), "N > 0");
else
  % default value from 'Design example' in article
  N = 4;
end

% Number of bits in the input data stream (including sign bit)
if exist('Bin', 'var')
  assertintegerscalar(Bin);
  assert((Bin > 0), "Bin > 0");
else
  % default value from 'Design example' in article
  Bin = 16;
end

% Number of bits in the output data stream (including sign bit)
if exist('Bout', 'var')
  assertintegerscalar(Bout);
  assert((Bout > 0), "Bout > 0");
else
  % default value from 'Design example' in article
  Bout = 16;
end

% Implementation of register widths reducing (at all stages):
%  true  - truncation
%  false - rounding
if exist('truncate_not_round', 'var')
  assert(isscalar(truncate_not_round) && islogical(truncate_not_round), ...
         "truncate_not_round is valid boolean");
else
  % default selection from 'Design example' in article
  truncate_not_round = true;
end

% First Bj values to override.
% Vector of length not exceeding total number of stages minus one
% (final last stage always reduces output according to Bout).
% Setting 'overrideB = []' with all other parameters defaults reproduces
% results from 'Design example' in article.
% Setting 'overrideB = zeros(1, 2*N)' skips register pruning completely.
if exist('overrideB', 'var')
  assert(   isempty(overrideB) ...
         || (   isvector(overrideB) ...
             && isnumeric(overrideB) ...
             && (length(overrideB) <= 2*N) ...
             && all(overrideB == round(overrideB)) ...
             && all((overrideB >= 0))), ...
         "overrideB is empty or valid non-negative integers vector of size not exceeding 2N" ...
        );
else
  % default is to force zero B₁ to avoid decreasing SNR even before filtering
  overrideB = 0;
end

clear assertintegerscalar;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate base parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Gmax = (R*M)^N;
Bmax = ceil(N*log2(R*M) + Bin - 1);
j_2np1 = 2*N+1; % total number of stages

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Allocate some vectors and initialize values to zero
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

h = cell(2*N, 1);
B = zeros(1, j_2np1);
E = zeros(1, j_2np1);
mu_j = zeros(1, j_2np1);
sigma_j = zeros(1, j_2np1);
D = zeros(1, j_2np1);
mu_tj = zeros(1, j_2np1);
F = zeros(1, j_2np1);
sigma_tj = zeros(1, j_2np1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate hj (impulse response coefficients of system
%               from jth stage up to and including the last stage)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for j = 1:2*N
  if (j <= N)
    hj_len = (R*M - 1)*N + j;
    h(j) = mat2cell(zeros(1, hj_len), 1);
    for k = 0:hj_len-1
      hjk = 0;
      for l = 0:floor(k/(R*M))
        hjk = hjk + ((-1)^l) * nchoosek(N, l) * nchoosek((N - j + k - R*M*l), (k - R*M*l));
      end
      h{j}(k+1) = hjk;
    end
  else
    hj_len = 2*N + 1 - j + 1;
    h(j) = mat2cell(zeros(1, hj_len), 1);
    for k = 0:hj_len-1
      h{j}(k+1) = ((-1)^k) * nchoosek((2*N + 1 - j), k);
    end
  end
end
clear hj_len hjk;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate "mean error gain" (simplified)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

D(1) = Gmax;
D(2:2*N) = 0;
D(j_2np1) = 1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate "variance error gain"
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for j = 1:2*N
  F(j) = sqrt(sum(h{j}.^2));
end
F(j_2np1) = 1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Set overriden values and last value of Bj
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

B(1:length(overrideB)) = overrideB;
B(j_2np1) = Bmax - Bout + 1;
B(j_2np1) = max([B(j_2np1) 0]); % output register width may exceed accumulation width

if (B(j_2np1) > 0) % no sense to prune registers when final output isn't reduced at all
  % Simplified pre-calculation of sigma_t parameter for last stage because
  % it needed to calculate Bj values for sections stages.
  sigma_tj(j_2np1) = sqrt((1/12)*((2^B(j_2np1))^2));
  % Calculate Bj values for stages in sections based on pre-calculated
  % parameters for last stage.
  for j = (length(overrideB)+1):(2*N) % skip overriden Bj values
    B(j) = floor(-log2(F(j)) + log2(sigma_tj(j_2np1)) + (1/2)*log2(6/N));
    B(j) = max([B(j) 0]); % there are cases when value may be negative
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate left parameters for all stages based on Bj values
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for j = 1:j_2np1
  if (B(j) == 0)
    E(j) = 0;
  else
    E(j) = 2^B(j);
  end
  if truncate_not_round
    mu_j(j) = (1/2)*E(j);
  else
    mu_j(j) = 0;
  end
  sigma_j(j) = sqrt((1/12)*(E(j)^2));
  mu_tj(j) = mu_j(j)*D(j);
  sigma_tj(j) = sqrt((sigma_j(j)^2)*(F(j)^2));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate final total parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

mu_t = mu_tj(1) + mu_tj(j_2np1); % simplified
mu_out = mu_t / (2^B(j_2np1));

sigma_t = sqrt(sum(sigma_tj.^2));
sigma_out = sigma_t / (2^B(j_2np1));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate extra parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% when R*M isn't power of two, output range will not occupy full Bout bits,
% i.e. max magnitude will be less than 2^(Bout-1)
max_mag_out = fix(Gmax * (2^(Bin-1)) / (2^B(j_2np1)));

error_rms_out = sqrt((mu_out^2) + (sigma_out^2));
error_db = 20*log10(error_rms_out / max_mag_out);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Prepare helper functions for results output
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if isoctave
  string = @(A, varargin) A;
  formattedDisplayText = @(X, varargin) disp(X);
  disptable = @(T) prettyprint(T);
else
  disptable = @(T) disp(T);
end
dispnamevalue = @(name,value) ...
    fprintf('%37s: %s\n', ...
            name, ...
            strtrim(formattedDisplayText(value, ...
                                         'UseTrueFalseForLogical', true) ...
                   ) ...
           );

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Display output parameters (properties) of filter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

disp('<strong>Input parameters:</strong>');
dispnamevalue('R', R);
dispnamevalue('M', M);
dispnamevalue('N', N);
dispnamevalue('Bin', Bin);
dispnamevalue('Bout', Bout);
dispnamevalue('Truncate (not round)', truncate_not_round);
if (isempty(overrideB))
  dispnamevalue('Override Bj', '-');
elseif (length(overrideB) == 1)
  dispnamevalue('Override B₁', overrideB);
else
  dispnamevalue(['Override B₁₋' sprintf('%c', (8320 + num2str(length(overrideB)) - '0'))], ...
                overrideB);
end
disp(' ');

disp('<strong>CIC decimator properties:</strong>');
disp(' ');

dispnamevalue('Num of bits growth due to gain', B(j_2np1));
dispnamevalue('Num of accum bits with no truncation', (Bmax+1));
dispnamevalue('Gmax', Gmax);
dispnamevalue('Bmax', Bmax);
dispnamevalue('μᴛ', sprintf('%i (%.5g out units)', mu_t, mu_out));
dispnamevalue('σᴛ', sprintf('%.5g (%.5g out units)', sigma_t, sigma_out));
dispnamevalue('Max magnitude, out units', max_mag_out);
dispnamevalue('Error RMS, out units', sprintf('%.5g (%.5g dB)', error_rms_out, error_db));
disp(' ');

disptable(table( ...
    (1:j_2np1)', B', (Bmax+1-B)', E', mu_j', D', mu_tj', sigma_j', F', sigma_tj', ...
    cellstr([ ...
        string(char('overriden'.*((1:2*N) <= length(overrideB))')); ...
        string(char('output reducing'*(B(j_2np1)>0))) ...
    ]), ...
    'VariableNames', ...
    cellstr([ ...
        "j (stage)"; "Bj"; "Accum width"; "Ej"; "μj"; "Dj"; "μᴛj"; "σj"; "Fj"; "σᴛj"; ...
        "comment" ...
    ]) ...
));
disp(' ');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Filter analysis
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

h_full = h{1};
if isoctave
  [~, signal_pkg_flag] = pkg("describe", "signal");
  if strcmpi(signal_pkg_flag{1}, 'Loaded')
    figure("name", "Frequency response");
    freqz(h_full, 1);
    figure("name", "Impulse response");
    impz(h_full, 1);
    figure("name", "Group/phase delay");
    grpdelay(h_full, 1);
    figure("name", "Poles and zeros");
    zplane(h_full, 1);
  end
  clear signal_pkg_flag;
else
  if exist('fvtool') %#ok<EXIST>
    fvtool(h_full, 1);
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Final clean-up
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear dispnamevalue;
if isoctave
  clear string formattedDisplayText;
end
clear disptable;

clear j k l isoctave;

It's an much improved version of script from code snippet at Computing CIC Filter Register Pruning Using Matlab [Updated].
I failed to register at DSPRelated.com to leave comments there, so I publish code here.