% 
% Aliasing Reduction in Clipped Signals
% Fabian Esqueda, Stefan Bilao and Vesa Välimäki
% Aalto University, Dept. of Signal Processing and Acoustics
% 
% MATLAB code used to implement Figs. 1 and 2
% 
% August 2015
% 
%%
clear all; close all; clc;

SR = 44100; %Sampling Rate
Ts = 1/SR;  %Sampling Period
f0 = 1245;  %Fundamental Frequency
 t = Ts*(0:SR);  %Time Vector
 L = 0.45;  %Clipping Threshold
OS = 10;   % Oversampling Rate
cycles = 2;  %Cycles to plot
 
% NONBANLIMITED CLIPPED SINE WAVE %
x = cos(2*pi*f0*t);
cl_sine = x;
cl_sine(x > L) = L;
cl_sine(x < -L) = -L;

% OVERSAMPLED VERSION (JUST FOR PLOT)
os_cl_sine = interp(x,OS);
os_sine = os_cl_sine;
os_cl_sine(os_cl_sine > L) = L;
os_cl_sine(os_cl_sine < -L) = -L;
% bl_sine = cl_bl_sine + (Cm(m) * cos(2*pi*f0*m*t));


% BANLIMITED CLIPPED SINE WAVE %
Ny = round(SR*0.5); % Nyquist Frequency
max_harm = 2*floor( ((SR/(2*f0))-1)/2 ) + 1; % Maximum harmonic to synthesize
% See "A Worst-Case Analysis for Analog-Quality (Alias-Free) Digital
% Dynamics Processing" by Dan Mapes-Riordan
tau = 0.5 - (asin(L)/pi);
cl_bl_sine = 0;
for m=1:ceil(max_harm/2)
    m = 2*m - 1;
    if m==1
        Cm(m) = 1 + 2*tau*( 2*L*sinc(m*tau) - sinc((m-1)*tau) - sinc((m+1)*tau) );
    else
        Cm(m) = 2*tau*( 2*L*sinc(m*tau) - sinc((m-1)*tau) - sinc((m+1)*tau) );
    end
    cl_bl_sine = cl_bl_sine + (Cm(m) * cos(2*pi*f0*m*t));
end


% PLOT BOTH SIGNALS AND THEIR SPECTRUMS
% Plot Parameters
figure(1)
set(1, 'Color', 'w');
fontsize = 16;
fontname = 'Times';
solidLineWidth = 2;
dottedLineWidth = 1;
ff = 0:(SR-1);
gap = OS;

% **************************************************** %
% NONBANDLIMITED WAVES
os_sine_samples = 1:length(os_sine);
subplot(2,2,1), plot(os_sine_samples(1:gap:end),os_sine(1:gap:end), 'k--','LineWidth',dottedLineWidth), hold on;
subplot(2,2,1), plot(os_cl_sine,'k','LineWidth',solidLineWidth);
set(gca,'fontsize',fontsize,'fontname',fontname);
xlabel({'Time (ms)','(a)'}');
set(gca,'XTick', [1, 0.5*OS*SR/f0, OS*SR/f0, 1.5*OS*SR/f0, 2*OS*SR/f0]); 
set(gca,'XTicklabel',[{0}, 0.4, 0.8, 1.2, 1.6]);
axis([1 cycles*OS*SR/f0 -1.2 1.2])
% **************************************************** %
w = chebwin(length(cl_sine),200);
cl_sine_freq = abs(fft((w.*cl_sine'),SR));
cl_sine_freq = db(cl_sine_freq/max(cl_sine_freq));
subplot(2,2,2), plot(ff,cl_sine_freq,'k','LineWidth',1);
hold on;
cl_bl_sine_freq = abs(fft(cl_bl_sine,SR));
cl_bl_sine_freq = db(cl_bl_sine_freq/max(cl_bl_sine_freq));
subplot(2,2,2), plot(ff,cl_bl_sine_freq,'ko','LineWidth',1);
set(gca,'fontsize',fontsize,'fontname',fontname);
xlabel({'Frequency (kHz)','(b)'}');
ylabel('Magnitude (dB)');
% set(gca,'XTick', [0, 5000, 10000, 15000, 20000]); 
% set(gca,'XTicklabel',{0, 5, 10, 15, 20});
% set(gca,'YTick', [-200, -160, -120, -80, -40, 0]); 
% set(gca,'YTicklabel',{-200, -160, -120, -80, -40, 0});
xlim([0 SR/2])
ylim([-160 10])
% **************************************************** %
% BANDLIMITED WAVES
os_sine_samples = 1:length(os_sine);
subplot(2,2,3), plot(os_sine_samples(1:gap:end),os_sine(1:gap:end), 'k--','LineWidth',dottedLineWidth), hold on;
subplot(2,2,3), plot(interp(cl_bl_sine,OS),'k','LineWidth',solidLineWidth);
set(gca,'fontsize',fontsize,'fontname',fontname);
xlabel({'Time (ms)','(c)'}');
set(gca,'XTick', [1, 0.5*OS*SR/f0, OS*SR/f0, 1.5*OS*SR/f0, 2*OS*SR/f0]); 
set(gca,'XTicklabel',[{0}, 0.4, 0.8, 1.2, 1.6]);
axis([1 cycles*OS*SR/f0 -1.2 1.2])
% **************************************************** %
cl_bl_sine_freq = abs(fft(cl_bl_sine,SR));
cl_bl_sine_freq = db(cl_bl_sine_freq/max(cl_bl_sine_freq));
subplot(2,2,4), stem(cl_bl_sine_freq,'k','LineWidth',1,'BaseValue',-200);
set(gca,'fontsize',fontsize,'fontname',fontname);
xlabel({'Frequency (kHz)','(d)'}');
ylabel('Magnitude (dB)');
set(gca,'XTick', [0, 5000, 10000, 15000, 20000]); 
set(gca,'XTicklabel',{0, 5, 10, 15, 20});
set(gca,'YTick', [-200, -160, -120, -80, -40, 0]); 
set(gca,'YTicklabel',{-200, -160, -120, -80, -40, 0});
xlim([0 SR/2])
ylim([-160 10])
% **************************************************** %

% print -deps2 Fig-trivial-vs-bandlimited

x = cos(2*pi*f0*t);
os_sine = interp(x,OS);

% PLOT BOTH SIGNALS AND THEIR SPECTRUMS
% Plot Parameters
figure(2)
set(2, 'Color', 'w');
fontsize = 16;
fontname = 'Times';
solidLineWidth = 2;
dottedLineWidth = 1;
ff = 0:(SR-1);
gap = OS;

os_sine_samples = 1:length(os_sine);
subplot(2,2,1), plot(0:gap:length(os_sine)-1,os_sine(1:gap:end), 'k--','LineWidth',dottedLineWidth), hold on;
subplot(2,2,1), plot(0:length(os_cl_sine)-1,os_cl_sine,'k','LineWidth',solidLineWidth);
subplot(2,2,1), plot(OS*(0:length(cl_sine)-1),cl_sine,'k.','LineWidth',solidLineWidth,'MarkerSize',16);
set(gca,'fontsize',fontsize,'fontname',fontname);
xlabel({'Time (ms)','(a)'}');
set(gca,'XTick', [1, OS*5, OS*10]);
set(gca,'XTicklabel',[{0}, 0.1, 0.2]);

set(gca,'YTick', [L-0.01 L L+0.01]); 

axis([0 OS*10 0.44 0.46])

hold off;

% os_sine = decimate(os_sine,OS);
% length(os_sine)
OS2 = 100;
% os_sine = interp(os_sine,OS2);

os_sine_samples = 1:length(os_sine);
subplot(2,2,2), plot(interp(x,OS2), 'k--','LineWidth',dottedLineWidth), hold on;
subplot(2,2,2), plot(interp(cl_bl_sine,OS2),'k','LineWidth',solidLineWidth);
subplot(2,2,2), plot(OS2*(0:length(cl_bl_sine)-1),cl_bl_sine,'k.','LineWidth',solidLineWidth,'MarkerSize',16);

set(gca,'fontsize',fontsize,'fontname',fontname);
xlabel({'Time (ms)','(b)'}');
set(gca,'XTick', [1, OS2*5, OS2*10]);
set(gca,'XTicklabel',[{0}, 0.1, 0.2]);
set(gca,'YTick', [L-0.01 L L+0.01]); 

axis([0 OS2*10 0.44 0.46])

% print -deps2 Fig-trivial-vs-bandlimited-detail
