dan@0: function [M,fCen,Mfull]=mel_matrix(fs,NbCh,Nfft,warp,fhigh,flow)
dan@0: % [M,fCen,Mfull]=mel_matrix(fs,NbCh,Nfft,warp,fhigh)
dan@0: % returns Nfft/2+1 - by - NbCh matrix of triangular weights
dan@0: % NbCh filterbanks linearly spaced on a MEL scale till fhigh
dan@0: % warp: default 1.0
dan@0: % fhigh: default=fs/2
dan@0: % Mfull: matrix without truncation to Nfft/2+1 points
dan@0: 
dan@0: if nargin<4,
dan@0:     warp=1;
dan@0: end
dan@0: if nargin<5,
dan@0:     fhigh=fs/2;
dan@0: end
dan@0: 
dan@0: LowMel=2595 * log10(1+flow/700);
dan@0: NyqMel=2595 * log10(1+fhigh/700);
dan@0: 
dan@0: StartMel=LowMel + (0:NbCh-1)/(NbCh+1)*(NyqMel-LowMel);
dan@0: fCen=warp*700*(10.^(StartMel/2595)-1);
dan@0: StartBin=round(Nfft/fs*fCen)+1;
dan@0: 
dan@0: EndMel=LowMel + (2:NbCh+1)/(NbCh+1)*(NyqMel-LowMel);
dan@0: EndBin=round(warp*Nfft/fs*700*(10.^(EndMel/2595)-1))+1;
dan@0: 
dan@0: TotLen=EndBin-StartBin+1;
dan@0: 
dan@0: LowLen=[StartBin(2:NbCh) EndBin(NbCh-1)]-StartBin+1;
dan@0: HiLen=TotLen-LowLen+1;
dan@0: 
dan@0: M=sparse([],[],[],ceil(warp*Nfft/2+1),NbCh,sum(TotLen));
dan@0: for k=1:NbCh,
dan@0:     M(StartBin(k):StartBin(k)+LowLen(k)-1,k)=(1:LowLen(k))'/LowLen(k);
dan@0:     M(EndBin(k)-HiLen(k)+1:EndBin(k),k)=(HiLen(k):-1:1)'/HiLen(k);
dan@0: end
dan@0: %TotWeight=sum(M,1);
dan@0: Mfull=M;
dan@0: M=M(1:Nfft/2+1,:);
dan@0: %WeightSum=full([sum(M,1);TotWeight]);