close all; clear all; clc;

% Accompanying code for the DAFx-16 Paper: Rounding Corners with BLAMP
% by F. Esqueda, V. Välimäki and S. Bilbao
%
% Code used to generate Fig. 5 - Antialiased Triangular Oscillator
%
% Code written by F. Esqueda
% Last modified: 06-09-2016
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %

% DISCLAIMER: This code is provided "as is", without warranty of any kind.

% Sampling rate
Fs = 44100;
Ts = 1/Fs; 

% Fundamental frequency
f0 = 1661;

% Phase increment
T0 = f0/Fs;

% Slope
m = 4*T0;

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

% Output vector
y = zeros(Fs,1);

% Control variables used to synthesize waveform in a loop 

yn = 0;
yn1 = 0;
yn2 = 0;
yn3 = 0;

flg = 0;
app_correction = 0;

% Initial phase value
phi = 0.25;

for n=1:Fs
    
    % Upward portions
    if phi<=0.5
        
        % Synthesize triangle
        yn = 4*phi - 1;
        
        % Trivial signal can be obtained from this variable
        x(n) = yn;
        
        % Detect corners
        if flg == -1
            app_correction = 1;
            flag = 1;
        elseif app_correction == 1;
            app_correction = 0;
            
            % Fractional delay required to center polyBLAMP
            d =  (phi-T0)/T0;

            % Compute polyBLAMP coefficients
            h0 = -d^5/120 + d^4/24 - d^3/12 + d^2/12 - d/24 + 1/120;
            h2 = -d^5/40 + d^4/24 + d^3/12 + d^2/12 + d/24 + 1/120;
            h1 = d^5/40 - d^4/12 + d^2/3 - d/2 + 7/30;
            h3 = d^5/120;

            % Apply correction, superimpose scaled polyBLAMP
            yn  =  yn + 2*m*h0;
            yn1 = yn1 + 2*m*h1;
            yn2 = yn2 + 2*m*h2;
            yn3 = yn3 + 2*m*h3;

        end
        
        flg = 1;
    
    % Downard portions
    elseif phi>0.5
        
        % Synthesize triangle
        yn = 3-4*phi;
        
        % Trivial signal can be obtained from this variable
        x(n) = yn;
        
        % Detect corners
        if flg == 1
            app_correction = 1;
        elseif app_correction == 1;
            app_correction = 0;
            
            % Fractional delay required to center polyBLAMP
            d =  (phi-0.5-T0)/T0;

            % Compute polyBLAMP coefficients
            h0 = -d^5/120 + d^4/24 - d^3/12 + d^2/12 - d/24 + 1/120;
            h2 = -d^5/40 + d^4/24 + d^3/12 + d^2/12 + d/24 + 1/120;
            h1 = d^5/40 - d^4/12 + d^2/3 - d/2 + 7/30;
            h3 = d^5/120;
            
            % Apply correction, superimpose scaled polyBLAMP
            yn  =  yn - 1*2*m*h0;
            yn1 = yn1 - 1*2*m*h1;
            yn2 = yn2 - 1*2*m*h2;
            yn3 = yn3 - 1*2*m*h3;
            
        end
        
        flg = -1;
    end
   
    % Output sample
    y(n) = yn3;
    
    % Update ALL control variables
    
    yn3 = yn2;
    yn2 = yn1;
    yn1 = yn;
    
    % Compute new phase value (modulo counter)
    phi = mod(phi+T0,1);
    
end

y = [y(4:end); 0; 0; 0];

%% PLOT
figure(1)
solidLineWidth = 2;
dottedLineWidth = 1;
markerSize = 2;
fontName = 'Times';
fontSize = 16;
w = chebwin(length(y),200)';

% Harmonics Vector
Ny = round(Fs*0.5); % Nyquist Frequency
max_harm = floor(Ny/f0); % Maximum harmonic to synthesize
harms = f0*(1:1:max_harm);

subplot(2,2,1)
plot(x,'k','LineWidth',dottedLineWidth,'MarkerSize',markerSize), hold on;
plot(x,'ok','LineWidth',solidLineWidth,'MarkerSize',markerSize,'MarkerFaceColor','k'), hold off;
set(gca,'fontsize',fontSize,'fontname',fontName);
axis([1 2*Fs/f0 -1.1 1.1])
set(gca,'XTick',[1 20 40 60])
set(gca,'XTickLabel',[{'0'}, '20', '40', '60'])
set(gca,'YTick',[-1 0 1])
xlabel({'Time (samples)','(a)'})

subplot(2,2,2)
XF = abs(fft(w.*x,Fs*1));
XF = db(XF/max(XF));
plot(0:Fs-1,XF,'k','LineWidth',dottedLineWidth), hold on;
plot(harms,XF(harms+1),'ok','LineWidth',dottedLineWidth), hold off;
axis([0 Fs/2 -100 5])
set(gca,'fontsize',fontSize,'fontname',fontName);
set(gca,'YTick',-100:20:0)
set(gca,'XTick',[0 5e3 10e3 15e3 20e3])
set(gca,'XTickLabel',[{'0'},'5', '10' '15' '20'])
xlabel({'Frequency (kHz)','(b)'})
ylabel('Magnitude (dB)')


subplot(2,2,3)
plot(y,'k','LineWidth',dottedLineWidth,'MarkerSize',markerSize), hold on;
plot(y,'ok','LineWidth',solidLineWidth,'MarkerSize',markerSize,'MarkerFaceColor','k'), hold off;
set(gca,'fontsize',fontSize,'fontname',fontName);
axis([1 2*Fs/f0 -1.1 1.1])
set(gca,'XTick',[1 20 40 60])
set(gca,'XTickLabel',[{'0'}, '20', '40', '60'])
set(gca,'YTick',[-1 0 1])
xlabel({'Time (samples)','(c)'})


subplot(2,2,4)
YF = abs(fft(w.*y',Fs*1));
YF = db(YF/max(YF));
plot(0:Fs-1,YF,'k','LineWidth',dottedLineWidth), hold on;
plot(harms,YF(harms+1),'ok','LineWidth',dottedLineWidth), hold off;
axis([0 Fs/2 -100 5])
set(gca,'fontsize',fontSize,'fontname',fontName);
set(gca,'YTick',-100:20:0)
set(gca,'XTick',[0 5e3 10e3 15e3 20e3])
set(gca,'XTickLabel',[{'0'},'5', '10' '15' '20'])
xlabel({'Frequency (kHz)','(d)'})
ylabel('Magnitude (dB)')

% EOF