function x = sample_DBS_signal(N,aDBS,fDBS,dDBS,fsmp)
  # Point-samples a symmetric square pulse signal
  # of frequency {fDBS} (in Hz), amplitude {aDBS},
  # and duty cycle {dDBS}, at {N} points,
  # with sampling frequency {fsmp} (in Hz).
  assert(fsmp > 2*fDBS) # Too complicated otherwise.  
  tDBS = 1/fDBS; # Period of DBS square wave.
  wDBS = dDBS*tDBS; # Width of each pulse
  rDBS = 0; # Time offset of current DBS period.
  tsmp = 1/fsmp; # Sampling interval.
  for ii = 1:N
    t = 0; # Time since start of sampling interval.
    if rDBS < wDBS
      v = + aDBS;
    elseif rDBS < tDBS/2
      v = 0;
    elseif rDBS < tDBS/2 + wDBS
      v = - aDBS;
    else
      v = 0;
    end
    x(ii,1) = v;
    rDBS += tsmp;
    while rDBS >= tDBS
      rDBS -= tDBS;
    end
  end
  x = x - sum(x)/N;