function [s32, s16, s8, s4] = welson(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;
  doq = f < fn;
  phaseInc = f16 / fs; % lowest freq
  invInc = 1.0 / phaseInc;
  invInc8 = 0.5 * invInc;
  invInc4 = 0.25 * invInc;
  invInc2 = 0.125 * invInc;
  invIncq = 0.0625 * invInc;
  pwq = max(min(1 - 50e-6 * f, 1), 0);
  phase = 0;
  
  q16 = zeros(1,fs);
  q8 = zeros(1,fs);
  q4 = zeros(1,fs);
  q2 = zeros(1,fs);
  q = zeros(1,fs);
  
  for i = 1 : fs
    phase = mod(phase + phaseInc, 1);
    phase2 = mod(phase - 0.5, 1);
    s1mp = 1.0 - phase;
    s1mp2 = 1.0 - phase2;
    q16(i) = 0.5 * sign(phase - 0.5) + 0.5 - blepDiff(phase * invInc) + blepDiff(phase2 * invInc) + blepDiff(s1mp * invInc) - blepDiff(s1mp2 * invInc);
    %q16(i) = 0.5 * sign(phase - 0.5) + 0.5;
    if (do8)
      p8 = mod(phase + phase, 1);
      phase2 = mod(p8 - 0.5, 1);
      s1mp = 1.0 - p8;
      s1mp2 = 1.0 - phase2;
      q8(i) = 0.5 * sign(p8 - 0.5) + 0.5 - blepDiff(p8 * invInc8) + blepDiff(phase2 * invInc8) + blepDiff(s1mp * invInc8) - blepDiff(s1mp2 * invInc8);
      %q8(i) = 0.5 * sign(p8 - 0.5) + 0.5;
    endif
    if (do4)
      p4 = mod(p8 + p8, 1);
      phase2 = mod(p4 - 0.5, 1);
      s1mp = 1.0 - p4;
      s1mp2 = 1.0 - phase2;
      q4(i) = 0.5 * sign(p4 - 0.5) + 0.5 - blepDiff(p4 * invInc4) + blepDiff(phase2 * invInc4) + blepDiff(s1mp * invInc4) - blepDiff(s1mp2 * invInc4);
      %q4(i) = 0.5 * sign(p4 - 0.5) + 0.5;
    endif
    if (do2)
      p2 = mod(p4 + p4, 1);
      phase2 = mod(p2 - 0.5, 1);
      s1mp = 1.0 - p2;
      s1mp2 = 1.0 - phase2;
      q2(i) = 0.5 * sign(p2 - 0.5) + 0.5 - blepDiff(p2 * invInc2) + blepDiff(phase2 * invInc2) + blepDiff(s1mp * invInc2) - blepDiff(s1mp2 * invInc2);
      %q2(i) = 0.5 * sign(p2 - 0.5) + 0.5;
    endif
    if (doq)
      pq = mod(p2 + p2, 1);
      phase2 = mod(pq - pwq, 1);
      s1mp = 1.0 - pq;
      s1mp2 = 1.0 - phase2;
      q(i) = 0.5 * sign(pq - pwq) + 0.5 - blepDiff(pq * invIncq) + blepDiff(phase2 * invIncq) + blepDiff(s1mp * invIncq) - blepDiff(s1mp2 * invIncq);
      %q(i) = 0.5 * sign(pq - pwq) + 0.5;
    endif
  endfor
  
  s32 = -(q16 + 0.454545454545455 * q8 + 0.256410256410256 * q4 + 0.121951219512195 * q2); % errore q8
  s16 = -(q8 + 0.454545454545455 * q4 + 0.256410256410256 * q2 + 0.178571428571429 * q);
  s8 = -(q4 + 0.454545454545455 * q2 + 0.370370370370370 * q);
  s4 = -(1.2 * q2 + 0.8 * q);
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