logfmap.m

function [M,N,G] = logfmap(I,L,H)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% logfmap.m
%
% [M,N] = logfmap(I,L,H)
%     Return a maxtrix for premultiplying spectrograms to map
%     the rows into a log frequency space.
%     Output map covers bins L to H of input
%     L must be larger than 1, since the lowest bin of the FFT
%     (corresponding to 0 Hz) cannot be represented on a 
%     log frequency axis.  Including bins close to 1 makes 
%     the number of output rows exponentially larger.
%     N returns the recovery matrix such that N*M is approximately I
%     (for dimensions L to H).
%     
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Convert base-1 indexing to base-0
L = L-1;
H = H-1;

ratio = (H-1)/H;
opr = round(log(L/H)/log(ratio));
ibin = L*exp([0:(opr-1)]*-log(ratio));

M = zeros(opr,I);

for i = 1:opr
  % Where do we sample this output bin?
  % Idea is to make them 1:1 at top, and progressively denser below
  % i.e. i = max -> bin = topbin, i = max-1 -> bin = topbin-1, 
  % but general form is bin = A exp (i/B)
  tt = pi*([0:(I-1)]-ibin(i));
  M(i,:) = (sin(tt)+eps)./(tt+eps);
end

% Normalize rows, but only if they are boosted by the operation
%% Fixup gain in bottom bins

G = ones(1,I);
G(1:(H+1)) = [0:H]./H;

% Inverse is just transpose plus scaling
N = (M.*repmat(G,opr,1))';