%% Tanh Clipping
%
% This script plots the mangitude response of a soft-clipped sinewave
% (tanh clipping) processed trivially and with the first- and second-order antialiasing
% forms proposed in the paper
% 
% "Antiderivative antialiasing for memoryless nonlinearities" by
% Stefan Bilbao, Fabián Esqueda, Julian D Parker and Vesa Välimäki
% IEEE Signal Processing Letters, Vol. 24, Issue 7, March 2017.
% 
% Second antiderivative of tanh is computed using lookup table. This function is not
% optimized for real-time implementation.
%
% Written by F. Esqueda. Last update: October 25, 2017.

close all; clear all; clc;

%% System Parameters

% Sampling Rate & Sampling Period
Fs = 44100;
Ts = 1/Fs;

% Time vector (1 second)
t = Ts*(0:Fs-1);

% FUNDAMENTAL FREQUENCY AND GAIN OF THE INPUT SIGNAL
'Fundamental Frequency:'
f0 = 16e3*rand(1) % Random fundamental between 0–16 kHz
A = 10;
x = A*cos(2*pi*f0*t);

% Output signals
y_first_order_aa = zeros(size(x));  % First-Order Antialiased
y_second_order_aa = zeros(size(x)); % Second-Order Antialiased

% Ill-conditioning threshold for methods
thresh = 10e-10;

%% Trivial Clipping

y_trivial = tanh(x);

%% First-Order Antialiasing

% State variables
xn1 = 0;
Fn1 = 0;

% Processing loop
for n=1:length(x)
    
    % Read input sample
    xn = x(n);
    
    % Compute first antiderivative
    Fn = log(cosh(xn));
    
    % Check for ill-conditioning
    if abs(xn-xn1)<thresh
        
        % Compute original function
        x_avg = 0.5*(xn + xn1);
        y_first_order_aa(n) = tanh(x_avg);
        
    else
        
        y_first_order_aa(n) = (Fn - Fn1) / (xn - xn1);
        
    end
    
    % Update states
    Fn1 = Fn;
    xn1 = xn;
    
end

%% Second-Order Antialiasing

% State variables
xn1 = 0;    xn2 = 0;
Fn1 = 0;
Gn1 = 0;    Gn2 = 0;

for n=1:length(x)

    % Read input sample
    xn = x(n);

    % Evaluate SECOND antiderivative
    Fn = Tanh_2nd_Antiderivative(xn);

    % First round of antialiasing
    if (abs(xn-xn1)<thresh)
        
        % First-order escape rule
        x_avg = 0.5*(xn + xn1);

        % Compute FIRST antiderivative
        Gn = log(cosh(x_avg));
        
    else
        
        % Antiderivative method
        Gn = (Fn - Fn1) / (xn - xn1);
        
    end

    if (abs(xn-xn2)<thresh)
                                
        % Second-order escape rule
        delta = 0.5 * (xn - 2*xn1 + xn2);
        if abs(delta) < thresh
            
            % Compute original function
            x_avg = 0.25*(xn + 2*xn1 + xn2);
            y_second_order_aa(n) = tanh(x_avg);
            
        else
            
            x_avg = 0.5*(xn + xn2);
                        
            % Compute FIRST antiderivative
            tmp1 = log(cosh(x_avg));
            
            % Compute SECOND antiderivative
            tmp2 = Tanh_2nd_Antiderivative(x_avg);
            
            % Output
            y_second_order_aa(n) = (2/delta)*(tmp1 + (Fn1 - tmp2)/delta);
           
        end

    else
        
        % Antiderivative Form
        y_second_order_aa(n) = 2 * (Gn - Gn1)/(xn - xn2);

    end

    % Update state variables
    Fn1 = Fn;
    Gn1 = Gn;
    xn2 = xn1;
    xn1 = xn;
    
end
%% PLOT RESULTS

fontSize = 12;
fontName = 'Times';

win = hann(length(x))';

YZ = abs(fft(win.*y_trivial,Fs));
max_ref = max(YZ);
YZ = db(YZ/max_ref);
subplot(311)
plot(0:Fs-1,YZ,'k')
set(gca,'FontSize',fontSize,'FontName',fontName);
set(gca,'XTick',0:5e3:20e3)
set(gca,'XTickLabel',[{'0'},'5','10','15','20']);
set(gca,'YTick',-100:20:20)
xlabel('Frequency (kHz)')
ylabel('Magnitude (dB)')
axis([0 Fs/2 -100 20])
title('Trivial')

YZ = abs(fft(win.*y_first_order_aa,Fs));
YZ = db(YZ/max_ref);
subplot(312)
plot(0:Fs-1,YZ,'k')
set(gca,'FontSize',fontSize,'FontName',fontName);
set(gca,'XTick',0:5e3:20e3)
set(gca,'XTickLabel',[{'0'},'5','10','15','20']);
set(gca,'YTick',-100:20:20)
xlabel('Frequency (kHz)')
ylabel('Magnitude (dB)')
axis([0 Fs/2 -100 20])
title('First-Order Antialiasing')

YZ = abs(fft(win.*y_second_order_aa,Fs));
YZ = db(YZ/max_ref);
subplot(313)
plot(0:Fs-1,YZ,'k')
set(gca,'FontSize',fontSize,'FontName',fontName);
set(gca,'XTick',0:5e3:20e3)
set(gca,'XTickLabel',[{'0'},'5','10','15','20']);
set(gca,'YTick',-100:20:20)
xlabel('Frequency (kHz)')
ylabel('Magnitude (dB)')
axis([0 Fs/2 -100 20])
title('Second-Order Antialiasing')

