function y = stairSaw(fs, f, N)
  fn = 0.5 * fs;
  if (f >= fn)
    y = zeros(1, fs);
    return;
  endif
  doHigh = N * f < fn;
  k = 1.0 / N;
  phaseInc = f / fs;
  invInc = 1.0 / phaseInc;
  invIncN = k * invInc;
  phase = 0;
  
  for i = 1 : fs
    phase = mod(phase + phaseInc, 1.0);
    s1mp = 1.0 - phase;
    y(i) = phase - blepDiff(phase * invInc) + blepDiff(s1mp * invInc);
    if (doHigh)
      pN = mod(N * phase, 1.0);
      s1mpN = 1.0 - pN;
      y(i) -= k * (pN - blepDiff(pN * invIncN) + blepDiff(s1mpN * invIncN));
    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