close all; clear; clc; format long;

% Implement four different diode clipper models and plot their
% spectrograms for a 10-second linear sine sweep from 20-5kHz.
% 
% Methods:
% 
% –Trivial (at audio rate)
% –Oversampled by 2
% –With First-order antiderivative antialising
% –With Second-order antiderivative antialiasing
% 
% Antiderivative antialiasing method is based on "Antiderivative 
% antialiasing for memoryless nonlinearities" by Bilbao et al., 2017.
% 
% Code written by F. Esqueda – 31 October 2017

Fs = 44.1e3;
Ts = 1/Fs;

 % Simulation Length
T = 10;

% Gain
A = 10;

% Time Vector
t = Ts*(0:T*Fs-1);

% Chirp (20-5kHz)
Vin = A*chirp(t,20,T/2,2500,'linear');

%% Trivial
Vout_Trivial = Diode_Clipper_Trivial(Vin);

%% AA_1
Vout_AA_1 = Diode_Clipper_AA_1(Vin);

%% AA_2
Vout_AA_2 = Diode_Clipper_AA_2(Vin);

%% Oversampled by 2
% New time vector and chirtp
t = (0.5*Ts)*(0:T*(2*Fs)-1);
Vin = A*chirp(t,20,T/2,2500,'linear');
Vout_OS = Diode_Clipper_Trivial(Vin);

%% PLOT

figure(1)
fontSize = 12;
fontName = 'Times';
set(groot, 'defaultAxesTickLabelInterpreter','latex');

%% Trivial

subplot(221)

% Compute STFT
[s,w,t] = spectrogram(Vout_Trivial,chebwin(1024),512,Fs,'yaxis');

% Normalize
s = s./max(max(abs(s)));    % Scale the maximum to be at 0 dB
s_small = (abs(s)<0.0001); % Find values smaller than -80 dB
s(s_small) = 0.0001;       % Force no value to be smaller than -80 dB

% Convert from angular to Hz
w = w*(0.5*Fs)/pi;          

% Adjust time
time = linspace(0,T,length(t));

% Plot spectrogram
imagesc(time,w,db(s))
axis xy
colormap(1-gray)
xlabel({'Time (s)','\quad\,\,(a)'},'Interpreter','latex')
set(gca,'FontSize',fontSize,'FontName',fontName);
set(gca,'YTick', 0:5e3:20e3)
set(gca,'YTickLabel', [{'0'},'5k','10k','15k','20k'])
ylabel('Frequency (Hz)')
ylim([0 20e3])
title('Audio Rate','Interpreter','latex','FontSize',11)

%% Oversampled by 2

subplot(222)

% Compute STFT
[s,w,t] = spectrogram(Vout_OS,chebwin(2048),1024,2*Fs,'yaxis');

% Normalize
s = s./max(max(abs(s)));    % Scale the maximum to be at 0 dB
s_small = (abs(s)<0.0001); % Find values smaller than -80 dB
s(s_small) = 0.0001;       % Force no value to be smaller than -80 dB

% Convert from angular to Hz
w = w*Fs/pi;          

% Adjust time
time = linspace(0,T,length(t));

% Plot spectrogram
imagesc(time,w,db(s))
axis xy
colormap(1-gray)
xlabel({'Time (s)','\quad\,\,(b)'},'Interpreter','latex')
set(gca,'FontSize',fontSize,'FontName',fontName);
set(gca,'YTick', 0:5e3:20e3)
set(gca,'YTickLabel', [{'0'},'5k','10k','15k','20k'])
ylabel('Frequency (Hz)')
ylim([0 20e3])
title('Oversampling by 2','Interpreter','latex','FontSize',11)

%% AA1

subplot(223)

% Compute STFT
[s,w,t] = spectrogram(Vout_AA_1,chebwin(1024),512,Fs,'yaxis');

% Normalize
s = s./max(max(abs(s)));    % Scale the maximum to be at 0 dB
s_small = (abs(s)<0.0001); % Find values smaller than -80 dB
s(s_small) = 0.0001;       % Force no value to be smaller than -80 dB

% Convert from angular to Hz
w = w*(0.5*Fs)/pi;          

% Adjust time
time = linspace(0,T,length(t));

% Plot spectrogram
imagesc(time,w,db(s))
axis xy
colormap(1-gray)
xlabel({'Time (s)','\quad\,\,(c)'},'Interpreter','latex')
set(gca,'FontSize',fontSize,'FontName',fontName);
set(gca,'YTick', 0:5e3:20e3)
set(gca,'YTickLabel', [{'0'},'5k','10k','15k','20k'])
ylabel({'Frequency (Hz)'})
ylim([0 20e3])
title('$1^{st}$-Ord. Antialiasing','Interpreter','latex','FontSize',11)


%% AA2

subplot(224)

% Compute STFT
[s,w,t] = spectrogram(Vout_AA_2,chebwin(1024),512,Fs,'yaxis');

% Normalize
s = s./max(max(abs(s)));    % Scale the maximum to be at 0 dB
s_small = (abs(s)<0.0001); % Find values smaller than -80 dB
s(s_small) = 0.0001;       % Force no value to be smaller than -80 dB

% Convert from angular to Hz
w = w*(0.5*Fs)/pi;          

% Adjust time
time = linspace(0,T,length(t));

% Plot spectrogram
imagesc(time,w,db(s))
axis xy
colormap(1-gray)
xlabel({'Time (s)','\quad\,\,(d)'},'Interpreter','latex')
set(gca,'FontSize',fontSize,'FontName',fontName);
set(gca,'YTick', 0:5e3:20e3)
set(gca,'YTickLabel', [{'0'},'5k','10k','15k','20k'])
ylabel({'Frequency (Hz)'})
ylim([0 20e3])
title('$2^{nd}$-Ord. Antialiasing','Interpreter','latex','FontSize',11)
