% ARBITRARY-ORDER IIR ANTIDERIVATIVE ANTIALIASING -
% P. P. La Pastina, S. D'Angelo, L. Gabrielli
% paper submitted to DAFx20in21
%
% Example code for the AA-IIR method using complex conjugate pole pairs and
% hard clipping. We use a sine sweep as input with gain > 1.
%
% Works with MATLAB only


Fs = 44100;
input_gain = 10;        

% TEST 1
stopbdB = 60;
stopbF = 0.61 * Fs; 
order = 10;
runsweep(Fs,input_gain,stopbdB,stopbF,order,'cheby2');

% TEST 2
stopbF = 0.45 * Fs;
order = 2;
runsweep(Fs,input_gain,'',stopbF,order,'butter');

% TRIVIAL METHOD
runsweep(Fs,input_gain,'','','','trivial');


%%
function runsweep(Fs,input_gain,stopbdB,cutoff,order,type)

    duration = 10*Fs;
    
    fv = linspace(20,Fs/2,duration) / Fs;
    n = 1:duration;
    sweep = input_gain * sin(pi.*fv.*n);
    x = sweep;
    
    if strcmp(type,'cheby2')
        Fcrads = 2*pi*cutoff / Fs;
        [z,p,k] = cheby2(order, stopbdB, Fcrads, 's');
        [b,a] = zp2tf(z,p,k);
        iirba = [a; b];
        [r,p,k] = residue(b,a);
        name = 'AA-IIR-2';
    elseif strcmp(type,'butter')
        Fcrads = 2*pi*cutoff / Fs;
        [z,p,k] = butter(order, Fcrads, 's');
        [b,a] = zp2tf(z,p,k);
        iirba = [a; b];
        [r,p,k] = residue(b,a);
        name = 'AA-IIR-1';
    elseif strcmp(type,'trivial');
        y = min(max(x, -1), 1);
        name = 'TRIVIAL';
    else
        error('unknown method');
    end
    
    if strcmp(type,'cheby2') || strcmp(type,'butter')
        y = 0*x;
        for o = 1:2:order % poles are in conjg pairs, must take only one for each
            ri = r(o);
            zi = p(o);
            y = y + AAIIR_hard_complex(x, ri, zi, 1e-9);    
        end
    end
    
    WSIZE = 1024;
    bh = blackmanharris(WSIZE/4);
    figure, spectrogram(y, bh, 16, 2*WSIZE, Fs, 'MinThreshold', -80, 'yaxis');
    colorMap = [linspace(1,0,256)', linspace(1,0,256)', linspace(1,0,256)'];
    colormap(colorMap); colorbar('off');
    pbaspect([2 1 1]);

    pbaspect([3 1 1]);
    annotation('textbox', [0.75, 0.5, 0.14, 0.05], 'String', name, 'BackgroundColor', 'white', 'HorizontalAlignment','center','VerticalAlignment','middle')

    if strcmp(type,'cheby2') || strcmp(type,'butter')
        [h,wout] = freqs(iirba(2,:),iirba(1,:),linspace(0.01, 2*2*pi, 2000));
        figure, plot(wout/(2*pi)*Fs /1000,mag2db(abs(h)),'k'); %/1000 for kHz tick notation
        xlim([0 wout(end)/(2*pi)*Fs /1000]);
        ylim([-90 10]);
        xlabel('Frequency [kHz]');
        ylabel('Magnitude [dB]');
        pbaspect([2 1 1]);
        line([Fs/2 /1000; Fs/2 /1000], [-100; 10],'Color','red');
        grid on;
        annotation('textbox', [0.75, 0.5, 0.14, 0.05], 'String', name, 'BackgroundColor', 'white', 'HorizontalAlignment','center','VerticalAlignment','middle')
    end

end