function [x,mn,mx]=melbankm(p,n,fs,fl,fh,w)
% 0002 %MELBANKM determine matrix for a mel-spaced filterbank [X,MN,MX]=(P,N,FS,FL,FH,W)
% 0003 %
% 0004 % Inputs:    p   number of filters in filterbank
% 0005 %        n   length of fft
% 0006 %        fs  sample rate in Hz
% 0007 %        fl  low end of the lowest filter as a fraction of fs (default = 0)
% 0008 %        fh  high end of highest filter as a fraction of fs (default = 0.5)
% 0009 %        w   any sensible combination of the following:
% 0010 %              't'  triangular shaped filters in mel domain (default)
% 0011 %              'n'  hanning shaped filters in mel domain
% 0012 %              'm'  hamming shaped filters in mel domain
% 0013 %
% 0014 %              'z'  highest and lowest filters taper down to zero (default)
% 0015 %              'y'  lowest filter remains at 1 down to 0 frequency and
% 0016 %               highest filter remains at 1 up to nyquist freqency
% 0017 %
% 0018 %               If 'ty' or 'ny' is specified, the total power in the fft is preserved.
% 0019 %
% 0020 % Outputs:    x     a sparse matrix containing the filterbank amplitudes
% 0021 %              If x is the only output argument then size(x)=[p,1+floor(n/2)]
% 0022 %              otherwise size(x)=[p,mx-mn+1]
% 0023 %        mn   the lowest fft bin with a non-zero coefficient
% 0024 %        mx   the highest fft bin with a non-zero coefficient
% 0025 %
% 0026 % Usage:    f=fft(s);            f=fft(s);
% 0027 %        x=melbankm(p,n,fs);        [x,na,nb]=melbankm(p,n,fs);
% 0028 %        n2=1+floor(n/2);        z=log(x*(f(na:nb)).*conj(f(na:nb)));
% 0029 %        z=log(x*abs(f(1:n2)).^2);
% 0030 %        c=dct(z); c(1)=[];
% 0031 %
% 0032 % To plot filterbanks e.g.    plot(melbankm(20,256,8000)')
% 0033 %
% 0034 
% 0035 
% 0036 %      Copyright (C) Mike Brookes 1997
% 0037 %      Version: $Id: melbankm.m,v 1.4 2007/05/04 07:01:39 dmb Exp $
% 0038 %
% 0039 %   VOICEBOX is a MATLAB toolbox for speech processing.
% 0040 %   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
% 0041 %
% 0042 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 0043 %   This program is free software; you can redistribute it and/or modify
% 0044 %   it under the terms of the GNU General Public License as published by
% 0045 %   the Free Software Foundation; either version 2 of the License, or
% 0046 %   (at your option) any later version.
% 0047 %
% 0048 %   This program is distributed in the hope that it will be useful,
% 0049 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
% 0050 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% 0051 %   GNU General Public License for more details.
% 0052 %
% 0053 %   You can obtain a copy of the GNU General Public License from
% 0054 %   http://www.gnu.org/copyleft/gpl.html or by writing to
% 0055 %   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
% 0056 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 if nargin < 6
   w='tz';
   if nargin < 5
     fh=0.5;
     if nargin < 4
       fl=0;
     end
   end
 end
 f0=700/fs;
 fn2=floor(n/2);
 lr=log((f0+fh)/(f0+fl))/(p+1);
 % convert to fft bin numbers with 0 for DC term
 bl=n*((f0+fl)*exp([0 1 p p+1]*lr)-f0);
 b2=ceil(bl(2));
 b3=floor(bl(3));
 if any(w=='y')
   pf=log((f0+(b2:b3)/n)/(f0+fl))/lr;
   fp=floor(pf);
   r=[ones(1,b2) fp fp+1 p*ones(1,fn2-b3)];
   c=[1:b3+1 b2+1:fn2+1];
   v=2*[0.5 ones(1,b2-1) 1-pf+fp pf-fp ones(1,fn2-b3-1) 0.5];
   mn=1;
   mx=fn2+1;
 else
   b1=floor(bl(1))+1;
   b4=min(fn2,ceil(bl(4)))-1;
   pf=log((f0+(b1:b4)/n)/(f0+fl))/lr;
   fp=floor(pf);
   pm=pf-fp;
   k2=b2-b1+1;
   k3=b3-b1+1;
   k4=b4-b1+1;
   r=[fp(k2:k4) 1+fp(1:k3)];
   c=[k2:k4 1:k3];
   v=2*[1-pm(k2:k4) pm(1:k3)];
   mn=b1+1;
   mx=b4+1;
 end
 if any(w=='n')
   v=1-cos(v*pi/2);
 elseif any(w=='m')
   v=1-0.92/1.08*cos(v*pi/2);
 end
 if nargout > 1
   x=sparse(r,c,v);
 else
   x=sparse(r,c+mn-1,v,p,1+fn2);
 end