function [s32, s16, s8, s4] = welsonApprox(fs, f)
  fn = 0.5 * fs;
  f16 = 0.0625 * f; % lowest freq
  if (f16 >= fs)
    s32 = zeros(1, fs);
    s16 = zeros(1, fs);
    s8 = zeros(1, fs);
    s4 = zeros(1, fs);
    return;
  endif
  do8 = 2.0 * f16 < fn;
  do4 = 4.0 * f16 < fn;
  do2 = 8.0 * f16 < fn;
  do1 = f < fn;
  phaseInc = f16 / fs;
  invInc = 1.0 / phaseInc;
  invInc8 = 0.5 * invInc;
  invInc4 = 0.25 * invInc;
  invInc2 = 0.125 * invInc;
  invInc1 = 0.0625 * invInc;
  pwq = max(min(1 - 50e-6 * f, 1), 0);
  p16 = 0;
  
  s32 = zeros(1,fs);
  s16 = zeros(1,fs);
  s8 = zeros(1,fs);
  s4 = zeros(1,fs);
  
  for i = 1 : fs
    p16 = mod(p16 + phaseInc, 1);
    s1mp = 1.0 - p16;
    s32(i) = -2.0 * (p16 - blepDiff(p16 * invInc) + blepDiff(s1mp * invInc));
    %s32(i) = -2.0 * p16;
    if (do8)
      p8 = mod(p16 + p16, 1);
      s1mp = 1.0 - p8;
      s16(i) = -2.0 * (p8 - blepDiff(p8 * invInc8) + blepDiff(s1mp * invInc8));
      %s16(i) = -2.0 * p8;
    endif
    if (do4)
      p4 = mod(p8 + p8, 1);
      s1mp = 1.0 - p4;
      s8(i) = -2.0 * (p4 - blepDiff(p4 * invInc4) + blepDiff(s1mp * invInc4));
      %s8(i) = -2.0 * p4;
    endif
    if (do2)
      p2 = mod(p4 + p4, 1);
      s1mp = 1.0 - p2;
      s4(i) = -2.0 * (p2 - blepDiff(p2 * invInc2) + blepDiff(s1mp * invInc2));
      %s4(i) = -2.0 * p2;
    endif
    if (do1)
      p1 = mod(p2 + p2, 1);
      s1mp = 1.0 - p1;
      phase2 = mod(p1 - pwq, 1);
      s1mp2 = 1.0 - phase2;
      d = blepDiff(s1mp * invInc1) - blepDiff(p1 * invInc1);
      %d = 0;
      x1 = p1 + d;
      x2 = x1 - 0.35 * sign(p1 - pwq) - 0.35 - 0.7 * (d + blepDiff(phase2 * invInc1) - blepDiff(s1mp2 * invInc1));
      %x2 = x1 - 0.35 * sign(p1 - pwq) - 0.35;
      s32(i) += 0.125 * x1;
      s16(i) += 0.25 * x2;
      s8(i) += 0.5 * x2;
      s4(i) += x2;
    endif
  endfor
endfunction

function y = blepDiff(d)
  y = 0 * d;
  y(d<2) = polyval([-0.04166666666666666, 0.3333333333333333, -0.9166666666666666, 1.0, -0.3333333333333333], d(d<2));
  y(d<1) = polyval([0.125, -0.3333333333333333, -0.25, 1.0, -0.5], d(d<1));
endfunction