function h = constq(winN, fs, bpo, loEdge, hiEdge)
  ## compute an approximation to the constant-Q transform with bpo bins
  ## per octave.  The fs and winN parameters are respectively the sample
  ## rate of the audio input and the FFT window size, while the optional
  ## loEdge and hiEdge parameters are cutoffs.  This routine produces a
  ## matrix as output which multiplies the linear FFT to produce the
  ## constant-Q transform approximation.

  df = fs/winN;
  N = fs/2;

  ## if loEdge and hiEdge are not provided, choose sensible values (low
  ## A minus a quarter-tone for the low cutoff, 8kHz or the Nyquist
  ## frequency, whichever is lower as the high cutoff.
  if(nargin < 4)
    loEdge = 53.43;
    hiEdge = min(N,8000);
  endif
  nbins = bpo * log2(hiEdge/loEdge);

  ## central frequencies of the logarithmic (constant-Q) bins
  cfs = loEdge * 2^(-1/(2*bpo)) * 2.^([1:nbins]/bpo);
  ## approximate the bandwidths as symmetrical for now
  bws = (cfs*2^(1/(2*bpo)) - cfs*2^(-1/(2*bpo)))/2;

  lofs = cfs - bws;
  hifs = cfs + bws;

  h = zeros(nbins,winN);
  for x = (1:winN/2)
    h(:,x+1) = ((df*x > lofs) & (df*x < hifs));
  endfor

end