www.pudn.com > asr.rar > MELCEPdwt_slow.M, change:2004-11-06,size:5713b
function c=melcepdwt(s,fs,w,nc,p,n,inc,fl,fh)
%MELCEPdwt Calculate the mel cepstrum of a signal C=(S,FS,W,NC,P,N,INC,FL,FH)
%DWT substitude the DCT
%
% Simple use: c=melcepst(s,fs) % calculate mel cepstrum with 12 coefs, 256 sample frames
% c=melcepst(s,fs,'e0dD') % include log energy, 0th cepstral coef, delta and delta-delta coefs
%
% Inputs:
% s speech signal
% fs sample rate in Hz (default 11025)
% nc number of dwt cepstral coefficients (default 15)
% n length of frame (default power of 2 <30 ms))
% p number of filters in filterbank (default floor(3*log(fs)) )
% inc frame increment (default n/2)
% fl low end of the lowest filter as a fraction of fs (default = 0)
% fh high end of highest filter as a fraction of fs (default = 0.5)
%
% w any sensible combination of the following:
%
% 'R' rectangular window in time domain
% 'N' Hanning window in time domain
% 'M' Hamming window in time domain (default)
%
% 't' triangular shaped filters in mel domain (default)
% 'n' hanning shaped filters in mel domain
% 'm' hamming shaped filters in mel domain
%
% 'p' filters act in the power domain
% 'a' filters act in the absolute magnitude domain (default)
%
% '0' include 0'th order cepstral coefficient
% 'e' include log energy
% 'd' include delta coefficients (dc/dt)
% 'D' include delta-delta coefficients (d^2c/dt^2)
%
% 'z' highest and lowest filters taper down to zero (default)
% 'y' lowest filter remains at 1 down to 0 frequency and
% highest filter remains at 1 up to nyquist freqency
%
% If 'ty' or 'ny' is specified, the total power in the fft is preserved.
%
% Outputs: c mel cepstrum output: one frame per row
%
% Copyright (C) Mike Brookes 1997
%
% Last modified Thu Jun 15 09:14:48 2000
%
% VOICEBOX is a MATLAB toolbox for speech processing. Home page is at
% http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
% Modified by LXL
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You can obtain a copy of the GNU General Public License from
% ftp://prep.ai.mit.edu/pub/gnu/COPYING-2.0 or by writing to
% Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin<2 fs=11025; end
if nargin<3 w='M'; end
if nargin<4 nc=15; end
if nargin<5 p=floor(3*log(fs)); end
if nargin<6 n=pow2(floor(log2(0.03*fs))); end
if nargin<9
fh=0.5;
if nargin<8
fl=0;
if nargin<7
inc=floor(n/2);
end
end
end
if any(w=='R')
z=enframe(s,n,inc);
elseif any (w=='N')
z=enframe(s,hanning(n),inc);
else
z=enframe(s,hamming(n),inc);
end
f=rfft(z.');
[m,a,b]=melbankm(p,n,fs,fl,fh,w);
pw=f(a:b,:).*conj(f(a:b,:));
pth=max(pw(:))*1E-6;
if any(w=='p')
y=log(max(m*pw,pth));
else
ath=sqrt(pth);
y=log(max(m*abs(f(a:b,:)),ath));
end
%c=rdct(y).';
jn=size(y,2);
c=[];
for j=1:jn
cc=[];
[c4,ln4]=wavedec(y(:,j),2,'bior3.5');
[c4t,in4t]=sort(abs(c4));
len4t=length(c4t);
cc=[cc,c4(in4t(len4t-7)),c4(in4t(len4t-6)),c4(in4t(len4t-5)),c4(in4t(len4t-4)),...
c4(in4t(len4t-3)),c4(in4t(len4t-2)),c4(in4t(len4t-1)),c4(in4t(len4t))];
[c8,ln8]=wavedec(y(:,j),3,'bior3.5');
f2ln8=ln8(1)+ln8(2);
[c8t,in8t]=sort(abs(c8(1:f2ln8)));
len8t=length(c8t);
cc=[cc,c8(in8t(len8t-3)),c8(in8t(len8t-2)),c8(in8t(len8t-1)),c8(in8t(len8t))];
[c16,ln16]=wavedec(y(:,j),4,'bior3.5');
f2ln16=ln16(1)+ln16(2);
[c16t,in16t]=sort(abs(c16(1:f2ln16)));
len16t=length(c16t);
cc=[cc,c16(in16t(len16t-1)),c16(in16t(len16t))];
[c32,ln32]=wavedec(y(:,j),5,'bior3.5');
f2ln32=ln32(1)+ln32(2);
[c32t,ln32t]=max(abs(c32(1:f2ln32)));
cc=[cc,c32(ln32t)];
c=[c;cc];
end
nf=size(c,1);
%nc=nc+1;
%if p>nc
% c(:,nc+1:end)=[];
%elseif p<nc
% c=[c zeros(nf,nc-p)];
%end
%if ~any(w=='0')
% c(:,1)=[];
% nc=nc-1;
%end
if any(w=='e')
% energy normalization
egy=sum(pw);
egymax=max(egy);
egy=egy/egymax;
c=[log(egy).' c];
% c=[log(sum(pw)).' c];
end
% calculate derivative
if any(w=='D')
vf=(4:-1:-4)/60;
af=(1:-1:-1)/2;
ww=ones(5,1);
cx=[c(ww,:); c; c(nf*ww,:)];
ncc=nc;
if any(w=='e')
ncc=nc+1;
end
vx=reshape(filter(vf,1,cx(:)),nf+10,ncc);
vx(1:8,:)=[];
ax=reshape(filter(af,1,vx(:)),nf+2,ncc);
ax(1:2,:)=[];
vx([1 nf+2],:)=[];
if any(w=='d')
c=[c vx ax];
else
c=[c ax];
end
elseif any(w=='d')
vf=(4:-1:-4)/60;
ww=ones(4,1);
cx=[c(ww,:); c; c(nf*ww,:)];
ncc=nc;
if any(w=='e')
ncc=nc+1;
end
vx=reshape(filter(vf,1,cx(:)),nf+8,ncc);
vx(1:8,:)=[];
c=[c vx];
end
if nargout<1
[nf,nc]=size(c);
t=((0:nf-1)*inc+(n-1)/2)/fs;
ci=(1:nc)-any(w=='0')-any(w=='e');
imh = imagesc(t,ci,c.');
axis('xy');
xlabel('Time (s)');
ylabel('Mel-cepstrum coefficient');
map = (0:63)'/63;
colormap([map map map]);
colorbar;
end