www.pudn.com > asr.rar > MELCEPdwtNorm_bior.M, change:2004-11-06,size:5874b
function c=melcepdwtnorm(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).';
nc=15;
jn=size(y,2);
c=[];
for j=1:jn
cc=[];
z=y(:,j)';
[c5,ln5]=wavedec(z,5,'bior3.5');
cd2=detcoef(c5,ln5,2);
[cd2t,in2t]=sort(abs(cd2));
len2=length(cd2t);
cc=[cc,cd2(in2t(len2)),cd2(in2t(len2-1)),cd2(in2t(len2-2)),cd2(in2t(len2-3)),...
cd2(in2t(len2-4)),cd2(in2t(len2-5)),cd2(in2t(len2-6)),cd2(in2t(len2-7))];
cd3=detcoef(c5,ln5,3);
[cd3t,in3t]=sort(abs(cd3));
len3=length(cd3t);
cc=[cc,cd3(in3t(len3)),cd3(in3t(len3-1)),cd3(in3t(len3-2)),cd3(in3t(len3-3))];
cd4=detcoef(c5,ln5,4);
[cd4t,in4t]=sort(abs(cd4));
len4=length(cd4t);
cc=[cc,cd4(in4t(len4)),cd4(in4t(len4-1))];
cd5=detcoef(c5,ln5,5);
[cd5t,in5t]=max(abs(cd5));
cc=[cc,cd5(in5t)];
c=[c;cc];
end
[nf,nc]=size(c);
%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
% c is normalizated
cMean=mean(c);
cStd=std(c,1);
mMean=zeros(nf,nc);
for j=1:nc
for i=1:nf
mMean(i,j)=cMean(j);
end
end
sStd=zeros(nf,nc);
for j=1:nc
for i=1:nf
sStd(i,j)=cStd(j);
end
end
c=(c-mMean)./sStd;
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
[nf,nc]=size(c);
% 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