% FUNCTION FBM = FBANK_MATRIX( P )
%
% by Guillaume Lathoud at IDIAP
%
% Returns a p.n_filters - by - p.win_points matrix
%
% Typical use of FBM:
%   buffered        = buffer( signal, 512, 256 );
%   preemphasized   = filter( [ 1 -0.97 ], 1, buffered );
%   windowed        = repmat( hamming( 512 ), 1, size( preemphasized,2)) .* preemphasized;
%   mag_spec        = abs( fft( windowed ) );
%   fb_coefficients = FBM * mag_spec;
% 
%
% Parameters:
% p.fs         ( default 16000 )
% p.win_points ( default 512 )
% p.f_min      ( default 300 )
% p.f_max      ( default 5000 )
% p.n_filters  ( default 40 )

function fbm = fbank_matrix( p )

if nargin < 1
  p = [];
end

p_default.fs = 16000;
p_default.win_points = 512;
p_default.f_min = 300;
p_default.f_max = 5000;
p_default.n_filters = 40;

p = fill_default( p, p_default );

if p.f_max > p.fs / 2
  error( 'f_max > fs / 2 !!!');
end

% Move to Mel scale
mel_min     = 2595 * log10( 1 + p.f_min  / 700 );
mel_max     = 2595 * log10( 1 + p.f_max  / 700 );

% Define trianguler filters: left, right and center
mel_bandwidth = 2 * ( mel_max - mel_min ) / ( 1 + p.n_filters );
mel_center    = mel_min + (1:p.n_filters) * (mel_bandwidth/2);
mel_left      = mel_center - mel_bandwidth / 2;
mel_right     = mel_center + mel_bandwidth / 2;

% Move back to linear frequency scale
f_center = 700 * ( 10 .^ ( mel_center / 2595 ) - 1 );
f_left   = 700 * ( 10 .^ ( mel_left   / 2595 ) - 1 );
f_right  = 700 * ( 10 .^ ( mel_right  / 2595 ) - 1 );


% Discretize
% f=0 <-> k=1
% f=fs <-> k=p.win_points
k_center = 1 + round( f_center / p.fs * ( p.win_points-1 ) );
k_left   = 1 + ceil( f_left   / p.fs * ( p.win_points-1 ) );
k_right  = 1 + floor( f_right  / p.fs * ( p.win_points-1 ) );

% Define filter for magnitude spectrum
n_left  = 1 + k_center - k_left;
n_right = k_right - k_center;

fbm = zeros( p.n_filters, p.win_points );
for a = 1:p.n_filters
  fbm( a, k_left(a):k_center(a) ) = (1:n_left(a)) / n_left(a);
  fbm( a, k_center(a)+1:k_right(a) ) = (n_right(a):-1:1) / (1+n_right(a));
end

% DEBUG
if nargout < 1
  figure; plot( [f_left; f_center; f_right], repmat([0 1 0].', 1,  p.n_filters) )
  k = max(k_right);
  hold on; plot( (0:k-1)/(p.win_points-1)*p.fs, fbm(:,1:k).', 'o--' );
  title( ['true filters (continuous) and discretized versions (dashed' ...
	  ' & circles)'] );
  
end

