% Bucked Brigade Device (BBD) chip simulation with BLEP for alias-suppression.
%
% This script simulates a N-stages BBD chip, controlled by a (possibly time-varying) clock.
% The script does not attempt to simulate charge leakage, nonlinearities and noise.
% Hp: the script assumes the input is bandlimited, no antialias filtering is done.
%
% Input arguments: 
% x: input signal
% Fs: sampling frequency
% N: stages of the BBD (remember that the delay will be N/2 samples at C_R)
% CR: clock rate, can be a vector (same length as x) or a scalar
%
% Output arguments:
% y: output signal
% bbd: ending state of the delay line
%
% Author: S. D'Angelo (Orastron Srl) with contributions from: L. Gabrielli (UNIVPM)
% (C) 2025

function [y, bbd] = bbd_blep(x, Fs, N, C_R)
  
  
  % control rate
  clockIncV = C_R ./ Fs;
  invIncV = 1.0 ./ clockIncV;
  
  % initial conditions
  clock = 0.0;
  bbd = zeros(1, N / 2);
  bbd_i = 0;
  x_z1 = 0;
  bbd_z1 = 0;
  bbd_z2 = 0;
  bbd_z3 = 0;
  
  % audio rate processing  
  y = zeros(1, length(x));
  for i = 1:length(x)
    if ~isscalar(C_R)
        clockInc = clockIncV(i);
        invInc = invIncV(i);
    else
        clockInc = clockIncV;
        invInc = invIncV;
    end
    bbd_in = x(i);
    
    clock = clock + clockInc;
    while (clock >= 1.0)
      bbd_z3 = bbd_z2;
      bbd_z2 = bbd_z1;
      bbd_z1 = bbd(bbd_i + 1);
      
      d = invInc * (1.0 - clock) + 1.0;
      bbd(bbd_i + 1) = d * (bbd_in - x_z1) + x_z1;
      
      bbd_i = mod(bbd_i + 1, N / 2);
      clock = clock - 1.0;
    end
    
    y(i) = bbd(bbd_i + 1);
    bbd_zm1 = bbd(mod(bbd_i + 1, N / 2) + 1);
    bbd_zm2 = bbd(mod(bbd_i + 2, N / 2) + 1);
    bbd_zm3 = bbd(mod(bbd_i + 3, N / 2) + 1);
    bbd_d1 = clock * invInc;
    bbd_d2 = bbd_d1 + invInc;
    bbd_d3 = bbd_d2 + invInc;
    bbd_dm1 = (1 - clock) * invInc;
    bbd_dm2 = bbd_dm1 + invInc;
    bbd_dm3 = bbd_dm2 + invInc;
    y(i) = y(i) - ( (bbd_z1 - y(i)) * blepDiff(bbd_d1) ...
      + (bbd_z2 - bbd_z1) * blepDiff(bbd_d2) ...
      + (bbd_z3 - bbd_z2) * blepDiff(bbd_d3) ...
      + (bbd_zm1 - y(i)) * blepDiff(bbd_dm1) ...
      + (bbd_zm2 - bbd_zm1) * blepDiff(bbd_dm2) ...
      + (bbd_zm3 - bbd_zm2) * blepDiff(bbd_dm3) );
    
    x_z1 = bbd_in;
    
    end
  
end

function y = blepDiff(d)
  y = 0 * d;
  y(d<2) = polyval([-0.04166666666666666, 0.3333333333333333, -0.9166666666666666, 1.0, -0.3333333333333333], d(d<2));
  y(d<1) = polyval([0.125, -0.3333333333333333, -0.25, 1.0, -0.5], d(d<1));
end
