function [y_trivial, y_poly, y_tri, y_spline, y_trig, y_FSM] = sine_hard_sync(f0, f1, fs, L)
  
##  list of parameters
##  
##  f0 = frequency (Hz) of the slave oscillator
##  f1 = frequency (Hz) of the master oscillator
##  fs = sampling frequency
##  L = duration (samples) of the output
##  
##  list of outputs
##  
##  y_trivial = no antialiasing
##  y_poly = antialiasing with polynomial kernel
##  y_tri = antialiasing with triangular kernel
##  y_spline = antialiasing with B-spline kernel
##  y_trig = antialiasing with trigonometric kernel
##  y_FSM = antialiasing with frequency shifting method
  
  % constants
  
  f0 = f0/fs; % we convert frequencies from Hz to 1/samples
  f1 = f1/fs;
  
  omega_0 = 2 * pi * f0;
  
  T = 1/f1; % period (samples) of the master oscillator
  
  % polynomial kernel data
  
  a_poly = [1, -1]; % polynomial coefficients, notice that they are shifted with respect to the paper notation
  N_poly = length(a_poly) - 1;
  epsilon_poly = 1;
  
  % polynomial kernel constants
  
  p1 = zeros(N_poly + 1, 1);
  p2 = zeros(N_poly + 1, 1);
  q1 = zeros(N_poly + 1, 1);
  q2 = zeros(N_poly + 1, 1);  
  
  for l = 0:N_poly
    for m = 0:floor((N_poly - l - 1)/2)
      p1(l + 1) += (-1)^(m + l + 1) * factorial_quotient(2 * m + l + 1, l) * a_poly(2 * m + l + 2)/omega_0^(2 * m + 2);
      q1(l + 1) += (-1)^m * factorial_quotient(2 * m + l + 1, l) * a_poly(2 * m + l + 2)/omega_0^(2 * m + 2);      
    endfor
    for m = 0:floor((N_poly - l)/2)
      p2(l + 1) += (-1)^(m + l) * factorial_quotient(2 * m + l, l) * a_poly(2 * m + l + 1)/omega_0^(2 * m + 1);
      q2(l + 1) += (-1)^m * factorial_quotient(2 * m + l, l) * a_poly(2 * m + l + 1)/omega_0^(2 * m + 1);
    endfor
  endfor
  
  p1 = flip(p1); % in order to use polyval, polynomial coefficients have to be sorted the other way round
  p2 = flip(p2);
  q1 = flip(q1);
  q2 = flip(q2);
    
  G = abs(polyval(q1, epsilon_poly) * cos(omega_0 * epsilon_poly) + polyval(q2, epsilon_poly) * sin(omega_0 * epsilon_poly) + polyval(p1, 0));
    
  % polynomial residual function
  
  function y = residual_poly(t)
    
    if (t >= -epsilon_poly && t <= 0)
      y = sin(omega_0 * T/2)/G * (cos(omega_0 * (t + T/2 + epsilon_poly)) * polyval(p1, -epsilon_poly) - sin(omega_0 * (t + T/2 + epsilon_poly)) * polyval(p2, -epsilon_poly) - cos(omega_0 * T/2) * polyval(p1, t) + sin(omega_0 * T/2) * polyval(p2, t));
    elseif (t > 0 && t <= epsilon_poly)
      y = sin(omega_0 * T/2)/G * (cos(omega_0 * (t + T/2 + epsilon_poly)) * polyval(p1, -epsilon_poly) - sin(omega_0 * (t + T/2 + epsilon_poly)) * polyval(p2, -epsilon_poly) - 2 * cos(omega_0 * (t + T/2)) * polyval(p1, 0) - cos(omega_0 * T/2) * polyval(q1, t) + sin(omega_0 * T/2) * polyval(q2, t)) + 2 * sin(omega_0 * T/2) * cos(omega_0 * (t + T/2));
    else
      y = 0;
    endif
    
  endfunction
  
  % triangular residual function
  
  function y = residual_tri(t)
    
    if (t >= -1 && t <= 0)
      y = sin(omega_0 * T/2)/(1 - cos(omega_0)) * (cos(omega_0 * (t + T/2 + 1)) + omega_0 * sin(omega_0 * T/2) * (1 + t) - cos(omega_0 * T/2));
    elseif (t > 0 && t <= 1)
      y = sin(omega_0 * T/2)/(1 - cos(omega_0)) * (cos(omega_0 * (t + T/2 + 1)) - 2 * cos(omega_0 * (t + T/2)) + omega_0 * sin(omega_0 * T/2) * (1 - t) + cos(omega_0 * T/2)) + 2 * sin(omega_0 * T/2) * cos(omega_0 * (t + T/2));
    else
      y = 0;
    endif
    
  endfunction
  
  % b-spline residual function

  function y = residual_spline(t)
    
    if (t >= -1.5 && t <= -0.5)
      y = sin(omega_0 * T/2)/(2^3.5 * (1 - cos(omega_0))^1.5) * (8 * sin(omega_0 * (t + (T + 3)/2)) + 4 * omega_0^2 * sin(omega_0 * T/2) * t^2 + (12 * omega_0^2 * sin(omega_0 * T/2) - 8 * omega_0 * cos(omega_0 * T/2)) * t + (9 * omega_0^2 - 8) * sin(omega_0 * T/2) - 12 * omega_0 * cos(omega_0 * T/2));
    elseif (t > -0.5 && t <= 0.5)
      y = sin(omega_0 * T/2)/(2^2.5 * (1 - cos(omega_0))^1.5) * ((4 * cos(omega_0 * (T + 3)/2) - 12 * cos(omega_0 * (T + 1)/2)) * sin(omega_0 * t) + (4 * sin(omega_0 * (T + 3)/2) - 12 * sin(omega_0 * (T + 1)/2)) * cos(omega_0 * t) - 4 * omega_0^2 * sin(omega_0 * T/2) * t^2 + 8 * omega_0 * cos(omega_0 * T/2) * t + (3 * omega_0^2 + 8) * sin(omega_0 * T/2));
    elseif (t > 0.5 && t <= 1.5)
      y = sin(omega_0 * T/2)/(2^3.5 * (1 - cos(omega_0))^1.5) * ((8 * cos(omega_0 * (T + 3)/2) - 24 * cos(omega_0 * (T + 1)/2) + 24 * cos(omega_0 * (T - 1)/2)) * sin(omega_0 * t) + (8 * sin(omega_0 * (T + 3)/2) - 24 * sin(omega_0 * (T + 1)/2) + 24 * sin(omega_0 * (T - 1)/2)) * cos(omega_0 * t) + 4 * omega_0^2 * sin(omega_0 * T/2) * t^2 - (12 * omega_0^2 * sin(omega_0 * T/2) + 8 * omega_0 * cos(omega_0 * T/2)) * t + (9 * omega_0^2 - 8) * sin(omega_0 * T/2) + 12 * omega_0 * cos(omega_0 * T/2));
    else
      y = 0;
    endif
    
    if (t >= 0 && t <= 1.5)
      y += 2 * sin(omega_0 * T/2) * cos(omega_0 * (t + T/2));
    endif
    
  endfunction
  
  % trigonometric kernel data
  
  a_trig = [0.355768, 0.487396, 0.144232, 0.012604]; % trigonometric coefficients, notice that they are shifted with respect to the paper notation. this is a Nuttall window.
  N_trig = length(a_trig) - 1;
  epsilon_trig = 1;
  
  % trigonometric kernel constants
  
  K = 0;
  for k = 0:N_trig
    K += (-1)^k * a_trig(k + 1)/(omega_0^2 * epsilon_trig^2 - k^2 * pi^2);
  endfor
  K *= omega_0 * epsilon_trig;
  
  % trigonometric residual function
  
  function y = residual_trig(t)
    
    if (abs(t) <= epsilon_trig)
      sum_cos = 0;
      sum_sin = 0;
      for k = 0:N_trig
        sum_cos += a_trig(k + 1) * cos(k * pi/epsilon_trig * t)/(omega_0^2 * epsilon_trig^2 - k^2 * pi^2);
        sum_sin += k * a_trig(k + 1) * sin(k * pi/epsilon_trig * t)/(omega_0^2 * epsilon_trig^2 - k^2 * pi^2);
      endfor
      y = sin(omega_0 * T/2)/(K * sin(omega_0 * epsilon_trig)) * (omega_0 * epsilon_trig * sin(omega_0 * T/2) * sum_cos + pi * cos(omega_0 * T/2) * sum_sin - K * sin(omega_0 * (t + epsilon_trig + T/2)));
    else
      y = 0;
    endif
    
    if (t >= 0 && t <= epsilon_trig)
      y += 2 * sin(omega_0 * T/2) * cos(omega_0 * (t + T/2));
    endif
    
  endfunction
    
  % FSM constants
  
  phi_1 = omega_0 * T;
  phi_2 = 0;
  
  omega_L = pi - omega_0;
  omega_H = pi + omega_0;
  
  % audio rate
  
  for j = 1:L
    
    j_red = mod(j, T);
    
    y_trivial(j) = sin(omega_0 * j_red);
    
    y_poly(j) = y_trivial(j);
    for k = -floor(epsilon_poly/T):1 + floor(epsilon_poly/T)
      y_poly(j) += residual_poly(j_red - k * T);
    endfor
    
    y_tri(j) = y_trivial(j);
    for k = -floor(1/T):1 + floor(1/T)
      y_tri(j) += residual_tri(j_red - k * T);
    endfor
    
    y_spline(j) = y_trivial(j);
    for k = -floor(1.5/T):1 + floor(1.5/T)
      y_spline(j) += residual_spline(j_red - k * T);
    endfor
    
    y_trig(j) = y_trivial(j);
    for k = -floor(epsilon_trig/T):1 + floor(epsilon_trig/T)
      y_trig(j) += residual_trig(j_red - k * T);
    endfor
    
    aj_red = abs(j_red);
    if (aj_red < 1e-9)
      delta_h = kaiser_cont(j_red, 4, 2) * 0.5 * abs(i/pi * log(-omega_H/omega_L));
    elseif (aj_red < 1)
      delta_h = kaiser_cont(j_red, 4, 2) * 0.5 * (sign(j_red) - i/pi * (exp_int(-omega_H * j_red) - exp_int(omega_L * j_red)));
    else
      delta_h = 0;
    endif
    delta_X = exp(i * omega_0 * j_red) * (exp(i * phi_2) - exp(i * phi_1));
    y_FSM(j) = y_trivial(j) - imag(delta_h * delta_X);
      
  endfor
  
endfunction


function y = factorial_quotient(k, l)
  
  % this function computes k!/l!, assuming 0 <= l <= k
  
  y = 1;
  for i = l + 1:k
    y *= i;
  endfor
  
endfunction


function y = kaiser_cont(x, alpha, L)
  
  % this function computes a kaiser window of parameter alpha and length L, in continuous time
  
  if (abs(x) <= 0.5 * L)
    y = besseli(0, pi * alpha * sqrt(1 - (2 * x/L)^2))/besseli(0, pi * alpha);
  else
    y = 0;
  endif
  
endfunction


function y = exp_int(omega)
  
  % this function computes the exponential integral Ei as defined in Zavalishin's paper, that adopts a different convention from Matlab
  
  y = cosint(abs(omega)) + i * sinint(omega);
  
endfunction
