function [Y,MX] = logfsgram(X, N, SR, WIN, NOV, FMIN, BPO)
% [Y,MX] = logfsgram(X, N, SR, WIN, NOV, FMIN, BPO)
%    Calculate a log-frequency spectrogram
%    X is input signal; N is parent FFT window; SR is the source samplerate.
%    WIN is actual window length within FFT, NOV is number of overlapping 
%    points between successive windows.
%    Optional FMIN is the lowest frequency to display (80Hz);
%    BPO is the number of bins per octave (12).
%    MX returns the nlogbin x nfftbin mapping matrix;
%    sqrt(MX'*(Y.^2)) is an approximation to the original FFT
%    spectrogram that Y is based on, suitably blurred by going 
%    through the log-F domain.
% 2004-03-30 dpwe@ee.columbia.edu $Header: /homes/dpwe/matlab/columbiafns/RCS/logfsgram.m,v 1.3 2004/04/01 22:39:02 dpwe Exp $

if nargin < 2
  N = 1024;
end
if nargin < 3
  SR = 8000;
end
if nargin < 4
  WIN = [];
end
if nargin < 5
  NOV = [];
end
if nargin < 6
  FMIN = 80;
end
if nargin < 7
  BPO = 12;
end

if isempty(WIN)
  WIN = N;
end
if isempty(NOV)
  NOV = WIN/2;
end

% Calculate underlying STFT
XX = specgram(X,N,SR,WIN,NOV);

% Construct mapping matrix

% Ratio between adjacent frequencies in log-f axis
fratio = 2^(1/BPO);

% How many bins in log-f axis
nbins = floor( log((SR/2)/FMIN) / log(fratio) );

% Freqs corresponding to each bin in FFT
fftfrqs = [0:(N/2)]*(SR/N);
nfftbins = N/2+1;

% Freqs corresponding to each bin in log F output
logffrqs = FMIN * exp(log(2)*[0:(nbins-1)]/BPO);

% Bandwidths of each bin in log F
logfbws = logffrqs * (fratio - 1);

% .. but bandwidth cannot be less than FFT binwidth
logfbws = max(logfbws, SR/N);

% Controls how much overlap there is between adjacent bands
ovfctr = 0.5475;   % Adjusted by hand to make sum(mx'*mx) close to 1.0

% Weighting matrix mapping energy in FFT bins to logF bins
% is a set of Gaussian profiles depending on the difference in 
% frequencies, scaled by the bandwidth of that bin
freqdiff = ( repmat(logffrqs',1,nfftbins) - repmat(fftfrqs,nbins,1) )./repmat(ovfctr*logfbws',1,nfftbins);
mx = exp( -0.5*freqdiff.^2 );
% Normalize rows by sqrt(E), so multiplying by mx' gets approx orig spec back
mx = mx ./ repmat(sqrt(2*sum(mx.^2,2)), 1, nfftbins);

% Perform mapping in magnitude-squared (energy) domain
y = sqrt( mx * (abs(XX).^2) );

% so, we lost phase information...

if nargout < 1
  imagesc([0 length(X)/SR],[1 nbins],20*log10(y));
  axis xy
  xlabel('Time');
  ylabel('Frequency');
  yt = get(gca,'YTick');
  for i = 1:length(yt)
    ytl{i} = sprintf('%.0f',logffrqs(yt(i)));
  end
  set(gca,'YTickLabel',ytl);
else
  Y = y;
  MX = mx;
end