function x = simulate_DBS_signal(N,aDBS,fDBS,dDBS,fsmp)
  # Creates a signal {x} that is a symmetric square pulse 
  # train of frequency 
  # {fDBS}, amplitude {aDBS}, and duty cycle {dDBS},
  # sampled at a frequency {fsmp} with triangular
  # sampling window.
  # Supersample two simple DBS signals:
  HK = 4
  K = 2*HK+1; # Supersampling factor.
  dDBS1 = dDBS;
  fDBS1 = fDBS;
  s1 = sample_DBS_signal(K*N,0.98,fDBS1,dDBS1,K*fsmp);
  figure 41; plot(s1(1:200,1)); 
  title(sprintf('raw simulated %.1f Hz %.3f duty', fDBS1, dDBS1));
  
  dDBS2 = 1.0*dDBS;
  fDBS2 = 2*fDBS;
  s2 = sample_DBS_signal(K*N,0.98,fDBS2,dDBS2,K*fsmp);
  figure 42; plot(s2(1:200,1)); 
  title(sprintf('raw simulated %.1f Hz %.3f duty', fDBS2, dDBS2));
  
  # Combine them alternating with frequency 20 Hz:
  sDBS1 = 0;         # Vertical shift of signal {s1}.
  sDBS2 = -0.35; # Vertical shift of signal {s2}
  M = K*(fsmp/20);
  s = zeros(K*N,1);
  for ii = 1:K*N
    kk = mod(ii-1+floor(2*M/3),M);
    if kk < 0.50*M
      s(ii,1) = s1(ii,1) + sDBS1;
    else
      s(ii,1) = 0.40*s2(ii,1) + sDBS2;
    end 
    s(ii,1) *= aDBS;
  end
  figure 43; plot(s(1:200,1)); 
  title('raw simulated combined')
  # Downsample with filter:
  x = zeros(N,1);
  for ii = 1:N
    jj = K*(ii-1) + HK + 1;
    wc = 1.0/(HK+1); # Weight of center tap.
    sum = wc*s(jj,1);
    for kk = 1:HK
      w = wc*(1 - kk/(HK+1));
      sum += w*(s(jj-kk,1) +s(jj+kk,1));
    end
    x(ii,1) = sum;
  end
  figure 44; plot(x(1:200,1)); 
  title('smoothed simulated combined')
