% ARBITRARY-ORDER IIR ANTIDERIVATIVE ANTIALIASING -
% P. P. La Pastina, S. D'Angelo, L. Gabrielli
% paper submitted to DAFx20in21
%
% AA-IIR method for the Hard clipper, simple complex conjugate pole pair
% B is the residue, beta is the pole, following the paper notation
% tol is a numerical tolerance to compute the simpler expression
%
% Works with MATLAB and GNU/Octave

function y = aaiir_hard_complex(x, B, beta, tol)

    % precalculating some values
    E = exp(beta);
    F = (E - 1)/beta;

    % initial conditions
    x_z1 = 0;
    Y_z1 = 0;
    y = 0 * x;

    % audio processing loop
    for n = 1:length(x)
        if (abs(x(n) - x_z1) < tol)
            I = hard_clipper(0.5 * (x(n) + x_z1)) * F;
        else
            I = int(x_z1, x(n), beta);
        end
        Y = E * Y_z1 + 2*B * I;
        y(n) = real(Y);

        x_z1 = x(n);
        Y_z1 = Y;
    end
end

% computing the integral
function z = int(x_z1, x, s)
    if (x_z1 < x)
        if (x <= -1)
            z = (1 - exp(s))/s;
        elseif (x > -1 && x <= 1 && x_z1 <= -1)
            z = ((x - x_z1) * exp(s * (x + 1)/(x - x_z1)) - (s + 1) * x + x_z1 - s * exp(s))/s^2;
        elseif (x > 1 && x_z1 <= -1)
            z = ((x - x_z1) * (exp(s * (x + 1)/(x - x_z1)) - exp(s * (x - 1)/(x - x_z1))) - s * (exp(s) + 1))/s^2;
        elseif (x_z1 > -1 && x <= 1)
            z = (exp(s) * (x + (s - 1) * x_z1) - (s + 1) * x + x_z1)/s^2;
        elseif (x_z1 > -1 && x_z1 <= 1 && x > 1)
            z = ((x_z1 - x) * exp(s * (x - 1)/(x - x_z1)) + exp(s) * (x + (s - 1) * x_z1) - s)/s^2;
        else
            z = (exp(s) - 1)/s;
        end
    else
        if (x_z1 <= -1)
            z = (1 - exp(s))/s;
        elseif (x_z1 > -1 && x_z1 <= 1 && x <= -1)
            z = ((x_z1 - x) * exp(s * (x + 1)/(x - x_z1)) + exp(s) * (x + (s - 1) * x_z1) + s)/s^2;
        elseif (x_z1 > 1 && x <= -1)
            z = ((x_z1 - x) * (exp((s * (x + 1))/(x - x_z1)) - exp((s * (x - 1))/(x - x_z1))) + s * (exp(s) + 1))/s^2;
        elseif (x > -1 && x_z1 <= 1)
            z = (exp(s) * (x + (s - 1) * x_z1) - (s + 1) * x + x_z1)/s^2;
        elseif (x > -1 && x <= 1 && x_z1 > 1)
            z = ((x - x_z1) * exp((s * (x - 1))/(x - x_z1)) - (s + 1) * x + x_z1 + s * exp(s))/s^2;
        else
            z = (exp(s) - 1)/s;
        end
    end
end

function y = hard_clipper(x)
    y = min(max(x, -1), 1);
end
