function S = spectrogram(x,ddt,dwd,fmin,fmax,fsmp)
  % Returns a grayscale image with the spectrogram of signal{x}.
  % Each column of the image is the power spectrum 
  % of a window of {dwd} consecutive samples. Successive columns
  % coorresponf to windows shifted by {ddt} samples.
  %
  % Assumes that the sampling frequency is {fsmp} (in Hz).  Each power spectrum
  % is truncated to the range {[fmin _ fmax]} (in Hz)
  %
  % Currently requires {dwd} to be even.
  assert(mod(dwd,2) == 0)
  assert(size(x,2) == 1)
  Nx = size(x,1);
  Nt = floor((Nx - dwd)/ddt) + 1 % Number of sliding window positions.
  dfnyq = dwd/2; % discrete Nyquist frequency. May be fractional.
  Np = floor(dfnyq) + 1; % Number of frequencies in window power spectrum.
  fstep = fsmp/dwd  % Frequency increment per row.
  
  % Define number of rows and max frequency of spectrogram
  dfmin = max([0, floor(fmin*dwd/fsmp)]) % Min discrete freq to show.
  fmin = dfmin*fsmp/dwd % Actual min frequency.
  
  dfmax = min([floor(dfnyq), ceil(fmax*dwd/fsmp)]) % Max discrete freq to show.
  fmax = dfmax*fsmp/dwd % Actual max frequency.
  
  Nf = dfmax-dfmin+1 % Number of discrete frequencies in spectrogram.
  
  S = zeros(Nf,Nt);
  u = cosine_window(dwd,fsmp);
  for kw = 1:Nt
    % Extract a windowful:
    kx = (kw - 1)*ddt + 1;
    w = x(kx:kx+dwd-1,1);
    W = fft(w .* u);
    WP = power_spectrum_from_fft(W);
    assert(size(WP,1) == Np)
    assert(size(WP,2) == 1)
    S(1:Nf,kw) = flipud(WP(dfmin+1:dfmax+1,1));
  end;
  