function y = diode_clipper_bilinear_AAIIR(x, fs, alpha)
    
  % x = input of the system
  % fs = sampling frequency
  % alpha = simple real pole of the AA-IIR filter
  % y = output with AA-IIR with compensation
  
  % constants
  
  % circuit components
  
  VT = 0.025;
  Is = 1e-15;
  R = 1e3;
  C = 33e-9;
  
  N = 5; % number of subintervals for numerical integration
  
  % alpha
  
  A = -alpha; % this residue is chosen in order to have unitary DC gain
  
  % coefficients of the direct filter, with the numerical approximation (Equation (13))

  S1 = 0;
  for k = 1:N - 1
    S1 += k/N * exp(alpha * (1 - k/N));
  endfor
  
  S2 = 0;
  for k = 1:N - 1
    S2 += (1 - k/N) * exp(alpha * (1 - k/N));
  endfor
  
  b0 = A * (S1 + 0.5)/N;
  b1 = A * (S2 + 0.5 * exp(alpha))/N;
  a1 = -exp(alpha);
  
  % coefficients of the inverse filter
  
  b0_inv = 1/b0;
  b1_inv = a1/b0;
  a1_inv = b1/b0;  
  
  % fs
  
  k1 = 1/(2 * fs * R * C + 1);
  k2 = Is * R/(2 * fs * R * C + 1);
  k3 = 1/(VT * (2 * fs * R * C + 1));
  k4 = (VT * (2 * fs * R * C + 1) * log((Is * R)/(VT * (2 * fs * R * C + 1))) + Is * R)/(VT * (2 * fs * R * C + 1));
   
  function y = f(x)
    y = k1 * x + sign(x) * (k2 - VT * omega(k3 * abs(x) + k4));
  endfunction

  % initial conditions

  xi_z1 = 0;
  y_z1 = 0;
  Dy_z1 = 0;
  u_z1 = 0;
  
  % audio rate
  
  for i = 1:length(x)
    
    xi = x(i) + R * C * (2 * fs * y_z1 + Dy_z1);
    xi_diff = xi - xi_z1;
    T0 = xi_diff/N;
    I = 0.5 * (f(xi) + f(xi_z1) * exp(alpha));
    for k = 1:N - 1
      I += f(xi_z1 + k * T0) * exp(alpha * (1 - k/N));
    endfor
    I /= N;
    u = exp(alpha) * u_z1 + A * I;
    y(i) = b0_inv * u + b1_inv * u_z1 - a1_inv * y_z1;

    xi_z1 = xi;
    Dy_z1 = 2 * fs * (y(i) - y_z1) - Dy_z1;
    y_z1 = y(i);
    u_z1 = u;
    
  endfor
      
endfunction


function y = omega(x)
  
  if (x <= 400)
    y = lambertw(exp(x));
  else
    y = x - log(x);
  endif
  
endfunction