function [Y,MX,F,T] = logfsgram(X, N, SR, WIN, NOV, FMIN, BPO)
% [Y,MX,F,T] = 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.
%    F returns the frequencies, in Hz, of each row in Y, and T
%    returns the time in secs of each column.
% 2004-03-30 dpwe@ee.columbia.edu $Header: $

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
  % No outputs, so plot (like specgram)
  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
  % Outputs requested, so no plot
  Y = y;
  MX = mx;
  F = logffrqs;
  T = [0:size(Y,2)-1]*(WIN-NOV)/SR;
end
