function [bu,eu] = cleanup_signal(a,fmin,fmax,fsmp,th2)
  % Cleans the signal {a} by filtering out 
  % narrow bands around peaks in the power spectrum,
  % and applying a bandpass filter with cutoffs {fmin,fmax}.
  % Returns the input signal minus the freq peaks {bu} and 
  % the error {eu}, and the band-filtered signal {cu
  % both apodized.
  close all
  assert(size(a,2)==1)
  N = size(a,1) % Number of samples.
  assert(mod(N,2) == 0)
  NP = floor(N/2) + 1 % Number of entries in power spectrum.
  
  % Find prominent peaks in the spectrum:
  fpmin = 19.0;  % Ignore peaks below this frequency.
  fpmax = 70.0;  % Ignore peaks above this frequency.
  % fps = find_peaks(a,th2,fpmin,fpmax,fsmp);
  fps = [[20.1465, 0.090, 1]; [29.3009, 0.020, 1]; [40.2935, 0.180, 1]];
  % fps = [[ 130.047, 2.663, 1.00]; [239.967, 4.541, 1.00]];
  % fps = [[ 130.047, 2.663, 1.00];];
  % fps = [[239.93, 4.541, 1.00];];
  NF = size(fps,1)
  for kk = 1:NF
    printf("%7.3f ± %7.3f\n", fps(kk,1), fps(kk,2))
  end
  
  % Apodize signal for cleaner spectrum:
  M = ceil(1*fsmp) % 1 second worth of samples.
  u = apodize_window(N,M,fsmp); % Smooth out first & last {M} samples.
  au = a .* u;
  AU = fft(au);
  AUP = power_spectrum_from_fft(AU);
  assert(size(AUP,1) == NP)
  
  % Remove spectrum peaks from full signal:
  bu = u .* cleanup_freqs(au,fps,fsmp);  
  
  % Bandpass-filter the original signal:
  fnyq = fsmp/2; % Nyquist freq.
  if (fmin > 0) || (fmax < fnyq)
    % Choose low and high frequencies of each shoulder:
    fmin_lo = 0.75*fmin
    fmin_hi = 1.25*fmin
    
    if fmax <= fnyq
      fmax_lo = 0.66667*fmax 
      fmax_hi = fmax
    else
      fmax_lo = fnyq + 1
      fmax_hi = fnyq + 1
    end
    W = bandpass_filter(N,fmin_lo,fmin_hi,fmax_lo,fmax_hi,fsmp);
    BU = fft(bu);
    BU = BU .* W;
    bu = u .* ifft(BU);
  end
  
  BU = fft(bu);
  BUP = power_spectrum_from_fft(BU);
  
  eu = au - bu;
  EU = fft(eu);
  EUP = power_spectrum_from_fft(EU);
  
  % figure  1; timeplot(au,1,N,fsmp,'');      title('raw full');
  % figure  2; timeplot(au,600,850,fsmp,'');  title('raw detail');
  figure  3; pwrplot(AUP,0,NP-1,fsmp);      title('raw spectrum');
  % figure  4; timeplot(bu,1,N,fsmp,'');      title('clean full');
  % figure  5; timeplot(bu,600,850,fsmp,'');  title('clean detail');
  figure  6; pwrplot(BUP,0,NP-1,fsmp);      title('clean spectrum');
  % figure  7; timeplot(eu,1,N,fsmp,'');      title('noise full');
  % figure  8; timeplot(eu,600,850,fsmp,'');  title('noise detail');
  figure  9; pwrplot(EUP,0,NP-1,fsmp);      title('noise spectrum');