## Copyright (C) 2014 Stefano D'Angelo
##
## Permission is hereby granted, free of charge, to any person obtaining a copy
## of this software and associated documentation files (the "Software"), to deal
## in the Software without restriction, including without limitation the rights
## to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
## copies of the Software, and to permit persons to whom the Software is
## furnished to do so, subject to the following conditions:
##
## The above copyright notice and this permission notice shall be included in
## all copies or substantial portions of the Software.
##
## THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
## IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
## FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
## AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
## LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
## OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
## THE SOFTWARE.

## -*- texinfo -*-
## @deftypefn {Function File} {@var{y} =} buchla_lp_gate (@var{fs}, @var{x}, @var{cv}, @var{offset}, @var{trim}, @var{fb}, @var{mode}, @var{fbnorm}, @var{nl})
##
## Simulate the Buchla lowpass-gate filter as described in
## @quotation
## J. Parker and S. D'Angelo, ``A Digital Model of the Buchla Lowpass-Gate,'' in
## @emph{Proc. 16th Intl. Conf. Digital Audio Effects (DAFx-13)}, pp. 278-285,
## Maynooth, Ireland, September 2013.
## @end quotation
##
## It produces the output vector @var{y} from the audio input vector @var{x}
## and the control voltage input @var{cv}.  @var{cv} can either be a scalar or a
## vector whose length is equal to the length of @var{x}.  @var{fs} is the
## sample rate (in Hz).
##
## @var{offset} and @var{trim} are, respectively, the control voltage offset and
## the trimmer (@math{R_6}) setting ratios, and they can either be scalars or
## vectors whose length is equal to the length of @var{x}.
##
## @var{mode} is a scalar or vector whose length is equal to the length of
## @var{x} and @var{cv} which can contain the values: 1 (``Lowpass'' mode),
## 2 (``VCA'' mode), or 3 (``Both'' mode).
##
## @var{nl} indicates whether to simulate feedback nonlinearity (true) or not
## (false).
##
## @var{fb} is the feedback gain, either as the normalized ratio
## (@math{a_{norm}}) if @var{fbnorm} is true and @var{mode} is 1, or otherwise
## as the absolute value (@math{a}). @var{fb} can either be a scalar or a vector
## whose length is equal to the length of @var{x}.
##
## @end deftypefn

## Author: Stefano D'Angelo <zanga.mail@gmail.com>
## Maintainer: Stefano D'Angelo <zanga.mail@gmail.com>
## Version: 1.0.0
## Keywords: buchla lowpass gate vca filter

function y = buchla_lp_gate(fs, x, cv, offset, trim, fb, mode, fbnorm, nl)

  ### Constant values

  # Components

  C1    = 1e-9;
  C2    = 220e-12;
  C3lp  = 4.7e-9;
  Cc    = 2e-9;
  Ralp  = 5e6;
  Ravca = 5e3;
  R1R2  = 10e3;
  R3    = 150e3;
  R4    = 470e3;
  R5    = 100e3;
  R6max = 20e3;
  R7    = 33e3;
  R8    = 4.7e3;
  R9    = 470;
  Vs    = 15;
  G     = 2e5;
  n     = 3.9696;
  VT    = 26e-3;
  VB    = 3.9;
  Ifmin = 10e-6;
  Ifmax = 40e-3;
  trise = 12e-3 / 5;
  tfall = 250e-3 / 5;

  # Other values

  A     = 3.464;
  B     = 1136.212;
  xi    = 600;
  k0    = 1.468e2;
  k1    = 4.920e-1;
  k2    = 4.167e-4;
  k3    = 7.391e-9;
  kL    = 6.386;
  gamma = 0.001;

  ### Input arguments

  # Checks

  if (nargin != 9)
    print_usage();
  elseif (!isscalar(fs) || !isnumeric(fs) || !isreal(fs) || fs <= 0)
    error("buchla_lp_gate: FS must be a positive real");
  elseif (isscalar(x) || !isvector(x) || !isnumeric(x) || !isreal(x))
    error("buchla_lp_gate: X must be a vector of real values");
  elseif (!isnumeric(cv) || !isreal(cv)
          || (!isscalar(cv) && !isvector(cv))
          || (!isscalar(cv) && isvector(cv)
              && length(cv) != length(x)))
    error(["buchla_lp_gate: CV must be a scalar or a vector of real values " ...
           "whose length is equal to the length of X"]);
  elseif (!isnumeric(offset) || !isreal(offset)
          || (!isscalar(offset) && !isvector(offset))
          || (!isscalar(offset) && isvector(offset)
              && length(offset) != length(x)))
    error(["buchla_lp_gate: OFFSET must be a scalar or a vector of real " ...
           "values whose length is equal to the length of X"]);
  elseif (!isnumeric(trim) || !isreal(trim)
          || (!isscalar(trim) && !isvector(trim))
          || (!isscalar(trim) && isvector(trim)
              && length(trim) != length(x)))
    error(["buchla_lp_gate: TRIM must be a scalar or a vector of real " ...
           "values whose length is equal to the length of X"]);
  elseif (!isnumeric(fb) || !isreal(fb)
          || (!isscalar(fb) && !isvector(fb))
          || (!isscalar(fb) && isvector(fb)
              && length(fb) != length(x)))
    error(["buchla_lp_gate: FB must be a scalar or a vector of real values " ...
           "whose length is equal to the length of X"]);
  elseif (!isvector(mode))
    error("buchla_lp_gate: MODE must be a scalar or vector");
  else
    for i = 1:length(mode)
      if (!isnumeric(mode) || (mode(i) != 1 && mode(i) != 2 && mode(i) != 3))
        error("buchla_lp_gate: MODE must contain 1, 2, or 3");
      endif
    endfor
  endif

  # Convert scalar parameters to vectors

  if (isscalar(cv))
    cv = repmat(cv, 1, length(x));
  endif
  if (isscalar(offset))
    offset = repmat(offset, 1, length(x));
  endif
  if (isscalar(trim))
    trim = repmat(trim, 1, length(x));
  endif
  if (isscalar(fb))
    fb = repmat(fb, 1, length(x));
  endif
  if (isscalar(mode))
    mode = repmat(mode, 1, length(x));
  endif

  # Enforce bounds

  cv     = limit(cv, -10, 10, "CV");
  offset = limit(offset, 0, 1, "OFFSET");
  trim   = limit(trim, 0, 1, "TRIM");
  if (!nl)
    if (fbnorm)
      fb(mode == 1) = limit(fb(mode == 1), -inf, 1, "FB(mode == 1)");
    else
      fb(mode == 1) = limit(fb(mode == 1), -inf, 1 + (2 * C1 + C2) / C3lp,
                            "FB(mode == 1)");
    endif
  endif

  ### Non-constant values

  # Components

  R2 = R1R2 * offset;
  R1 = R1R2 - R2;
  R6 = R6max * trim;

  C3 = zeros(1, length(x));
  C3(mode == 1) = C3lp;
  Ra = repmat(Ralp, 1, length(x));
  Ra(mode == 2) = Ravca;

  # Other values

  alpha = 1 + (R6 + R7) * (1 / R3 + 1 / R5);
  beta  = (1 ./ alpha - 1) ./ (R6 + R7) - 1 / R8;

  ### Control path

  # High-shelving (see Errata)

  kVb1 = (1 + 2 * Cc * fs * (R4 + R5)) / (1 + 2 * Cc * fs * R4);
  kVb2 = (1 - 2 * Cc * fs * (R4 + R5)) / (1 + 2 * Cc * fs * R4);
  kVb3 = (1 - 2 * Cc * fs * R4) / (1 + 2 * Cc * fs * R4);

  Vb = filter([kVb1 kVb2], [1 kVb3], cv);
  Ia = Vb / R5 + Vs ./ (R3 * (1 + R1 ./ R2));

  # Static nonlinearity

  kIa      = xi * alpha * n * VT ./ (G * (R6 + R7 - 1 ./ (alpha .* beta)));
  c1       = Ia <= -kIa;
  c3       = Ia >= kIa;
  c2       = !(c1 | c3);
  V3       = zeros(1, length(x));
  V3(c1)   = -Ia(c1) ./ (alpha(c1) .* beta(c1));
  V3(c3)   = kL * alpha(c3) / G * n * VT - Ia(c3) .* (R6(c3) + R7);
  warg     = zeros(1, length(x));
  w        = zeros(1,length(x));
  warg(c2) = G * Ia(c2) .* (R6(c2) + R7 - 1 ./ (alpha(c2) .* beta(c2))) ...
             ./ (alpha(c2) * n * VT);
  w(c2)    = k0 + warg(c2) .* (k1 + warg(c2) .* (k2 + warg(c2) * k3));
  V3(c2)   = -alpha(c2) / G * n * VT .* w(c2) ...
             - Ia(c2) ./ (alpha(c2) .* beta(c2));

  kIa2   = alpha .* (Ifmin - beta .* V3);
  kIa3   = VB ./ (R6 + R7);
  kIa4   = (gamma * G * VB + alpha * R9 .* (VB * beta + Ifmax)) ...
           ./ (gamma * G * (R6 + R7) + R9);
  c1     = Ia <= kIa2;
  c2     = Ia <= kIa3 & !c1;
  c4     = Ia >= kIa4;
  c3     = !(c4 | c1 | c2);
  If     = zeros(1, length(x));
  If(c1) = Ifmin;
  If(c2) = beta(c2) .* V3(c2) + Ia(c2) ./ alpha(c2);
  If(c3) = gamma * G * (Ia(c3) .* (R6(c3) + R7) - VB) ./ (alpha(c3) * R9) ...
           - beta(c3) * VB + Ia(c3) ./ alpha(c3);
  If(c4) = Ifmax;

  Rf = A ./ (If .^ 1.4) + B;

  ### Vactrol model

  Rfz1       = [0 Rf(1:length(x) - 1)];
  dRf        = Rf - Rfz1;
  t          = repmat(tfall, 1, length(x));
  t(dRf < 0) = trise;
  k1         = 1 ./ (1 + 2 * fs * t);
  k2         = k1 .* (2 * fs * t - 1);

  Rfv   = zeros(1, length(x));
  Rfvz1 = 0;
  for i = 1:length(x)
    Rfv(i) = k1(i) * (Rf(i) + Rfz1(i)) + k2(i) * Rfvz1;
    Rfvz1  = Rfv(i);
  endfor

  ### Audio path

  # Denormalize feedback gain

  a = fb;
  if (fbnorm)
    a(mode == 1) .*= (2 * C1 * Ra(mode == 1) ...
                      + (C2 + C3(mode == 1)) ...
                      .* (Ra(mode == 1) + Rfv(mode == 1))) ...
                     ./ (C3(mode == 1) .* Ra(mode == 1));
  endif

  # Filter

  Rfinv = 1 ./ Rfv;
  a1    = Rfinv / C1;
  a2    = -(Rfinv + 1 ./ Ra) / C1;
  b1    = Rfinv / C2;
  b2    = -2 * b1;
  b3    = b1;
  b4    = C3 / C2;
  d1    = a;
  d2    = -1;

  koo   = 1 ./ (1 - a2 / (2 * fs));
  kox   = a1 / (2 * fs) .* koo;
  kdo   = 2 * fs * d1;
  kdx   = 2 * fs * d2;
  ksdo  = -4 * fs * d1;
  ksdx  = -4 * fs * d2;
  ksxi  = b1 / fs;
  ksxx  = b2 / fs;
  ksxo  = b3 / fs;
  ksxd  = b4 / fs;
  ksox  = a1 / fs;
  ksoo  = a2 / fs;

  sd   = 0;
  sx   = 0;
  so   = 0;
  y    = zeros(1, length(x));
  if (nl)
    yoz1 = 0;
    for i = 1:length(x)
      t    = tanh(yoz1);
      t2m1 = t * t - 1;

      # see Errata
      yxden  = a2(i) * b2(i) - a1(i) * b3(i) ...
               - 2 * fs * (a2(i) + b2(i) ...
                           + b4(i) * (a1(i) * d1(i) * t2m1 - a2(i) * d2)) ...
               - 4 * fs * fs * (b4(i) * d2 - 1);
      yxnom  = 2 * fs - a2(i);
      kyx    = yxnom / yxden;
      kxo    = 2 * fs * (-2 * fs * b4(i) * d1(i) * t2m1 + b3(i)) / yxden;
      kxi    = kyx * b1(i);
      kxd    = kyx * b4(i);
      kxx    = kyx * 2 * fs;
      kxt    = kxx * b4(i) * d1(i);
      kxyoz1 = kxt * t2m1;

      yx = kxo * so + kxi * x(i) + kxd * sd + kxx * sx + kxt * t ...
           + kxyoz1 * yoz1;
      yo   = koo(i) * so + kox(i) * yx;
      yonl = t - (yo - yoz1) * t2m1;
      yd   = sd + kdo(i) * yonl + kdx * yx;

      sd   = -sd + ksdo(i) * yonl + ksdx * yx;
      sx   = sx + ksxi(i) * x(i) + ksxx(i) * yx + ksxo(i) * yo + ksxd(i) * yd;
      so   = so + ksox(i) * yx + ksoo(i) * yo;
      yoz1 = yo;

      y(i) = yo;
    endfor
  else
    # see Errata
    yxden = a2 .* b2 - a1 .* b3 ...
            - 2 * fs * (a2 + b2 + b4 .* (a1 .* d1 - a2 * d2)) ...
            - 4 * fs * fs * (b4 * d2 - 1);
    yxnom = 2 * fs - a2;
    kyx   = yxnom ./ yxden;
    kxo   = 2 * fs * (2 * fs * b4 .* d1 + b3) ./ yxden;
    kxi   = kyx .* b1;
    kxd   = kyx .* b4;
    kxx   = kyx * 2 * fs;

    for i = 1:length(x)
      yx = kxo(i) * so + kxi(i) * x(i) + kxd(i) * sd + kxx(i) * sx;
      yo = koo(i) * so + kox(i) * yx;
      yd = sd + kdo(i) * yo + kdx * yx;

      sd = -sd + ksdo(i) * yo + ksdx * yx;
      sx = sx + ksxi(i) * x(i) + ksxx(i) * yx + ksxo(i) * yo + ksxd(i) * yd;
      so = so + ksox(i) * yx + ksoo(i) * yo;

      y(i) = yo;
    endfor
  endif

endfunction

function y = limit(x, low, up, name)
  y = x;
  if (min(y) < low || max(y) > up)
    warning(["buchla_lp_gate: limiting " name " to [" num2str(low) ", " ...
             num2str(up) "]"]);
    y(y < low) = low;
    y(y > up)  = up;
  endif
endfunction
