ACTIVLEV.M
上传用户:sy_best
上传日期:2013-01-11
资源大小:275k
文件大小:6k
源码类别:

语音合成与识别

开发平台:

Matlab

  1. function [lev,fso]=activlev(sp,fs,mode)
  2. %ACTIVLEV Measure active speech level as in ITU-T P.56 [LEV,FSO]=(SP,FS,MODE)
  3. %
  4. %Inputs: SP is the speech signal (with better than 20dB SNR)
  5. %        FS is the sample frequency in Hz (see also FSO below)
  6. %        MODE is any combination of the following:
  7. %            r - raw: do not apply the input filter 0.2 to 5.5 kHz
  8. %            a - all: use all the samples rather than subsampling to 694Hz
  9. %            d - give outputs in dB rather than power
  10. %            l - give the "long term level" as well as the active speech level
  11. %Outputs:
  12. %        LEV gives the speech level in units of power (or dB if mode='d')
  13. %            if mode='l' is specified, LEV is a row vector with the "long term
  14. %            level" as its second element (this is just the mean power)
  15. %        FSO is a column vector of intermediate information that allows
  16. %            you to process a speech signal in chunks. Thus:
  17. %
  18. %            fso=fs; for i=1:inc:nsamp, [lev,fso]=activlev(sp(i:i+inc-1),fso,mode); end
  19. %
  20. %            is equivalent to:                lev=activlev(sp(1:nsamp),fs,mode)
  21. %
  22. %            but is much slower. The two methods will not give identical results
  23. %            because they will use slightly different threshods.
  24. %For completeness we list here the contents of the FSO array:
  25. %
  26. %     1 : sample frequency
  27. %     2 : subsampling increment
  28. %     3 : hangover increment
  29. %   4-6 : smoothing filter coefs
  30. %  7-12 : 200Hz HP filter numerator
  31. % 13-18 : 200Hz HP filter denominator
  32. % 19-24 : 5.5kHz LP filter numerator
  33. % 25-30 : 5.5kHz LP filter denominator
  34. % 31-33 : count,sum and sum of squares of filter speech
  35. %    34 : upper limit of largest amplitude bin
  36. %    35 : samples to skip for subsampling
  37. % 36-37 : smoothing filter state
  38. % 38-42 : 200Hz HP filter state
  39. % 43-47 : 5.5kHz LP filter state
  40. % 48-55 : thresholded sample counts
  41. % 56-63 : hangover counts
  42. %      Copyright (C) Mike Brookes 1998
  43. %
  44. %      Last modified Thu Apr 30 17:28:29 1998
  45. %
  46. %   VOICEBOX home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
  47. %
  48. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  49. %   This program is free software; you can redistribute it and/or modify
  50. %   it under the terms of the GNU General Public License as published by
  51. %   the Free Software Foundation; either version 2 of the License, or
  52. %   (at your option) any later version.
  53. %
  54. %   This program is distributed in the hope that it will be useful,
  55. %   but WITHOUT ANY WARRANTY; without even the implied warranty of
  56. %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  57. %   GNU General Public License for more details.
  58. %
  59. %   You can obtain a copy of the GNU General Public License from
  60. %   ftp://prep.ai.mit.edu/pub/gnu/COPYING-2.0 or by writing to
  61. %   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
  62. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  63. if nargin<3 mode='0'; else mode = [mode '0']; end
  64. lev=zeros(1,1+any(mode=='l'));
  65. thf=0.5;
  66. nb=8;
  67. fso=[fs;zeros(64-1-length(fs),1)];
  68. if length(fs)==1
  69.    ni=max(floor(fs*all(mode~='a')/694),1);
  70.    ti=ni/fs;
  71.    nh=ceil(0.2/ti)+1;
  72.    g=exp(-ti/0.03);
  73.    % s-plane zeros and poles of high pass 5'th order chebychev2 filter with -0.25dB at w=1
  74.    szp=[0.37843443673309i 0.23388534441447i; -0.20640255179496+0.73942185906851i -0.54036889596392+0.45698784092898i]; 
  75.    szp=[[0; -0.66793268833792] szp conj(szp)];
  76.    zl=2./(1-szp*tan(200*pi/fs))-1;
  77.    al=real(poly(zl(2,:)));
  78.    bl=real(poly(zl(1,:)));
  79.    sw=1-2*rem(0:5,2).';
  80.    bl=bl*(al*sw)/(bl*sw);
  81.    fso(2:19-1)=[ni;nh;[1; -2*g; g^2]/(1-2*g+g^2);bl(:);al(:)];
  82.    if fs>14000
  83.       zh=2./(szp/tan(5500*pi/fs)-1)+1;
  84.       ah=real(poly(zh(2,:)));
  85.       bh=real(poly(zh(1,:)));
  86.       bh=bh*sum(ah)/sum(bh);
  87.       fso(19:19+11)=[bh(:);ah(:)];
  88.    end
  89. end
  90. % need to calculate offset correctly
  91. % subsample the signal
  92. ns=length(sp);
  93. if ns
  94.    % input filter goes here
  95.    if all(mode~='r')
  96.       [sp,fso(38:42)]=filter(fso(7:12),fso(13:18),sp(:),fso(38:42));
  97.       if fso(25)
  98.          [sp,fso(43:47)]=filter(fso(19:24),fso(25:30),sp,fso(43:47));
  99.       end
  100.    end
  101.    sp=abs(sp(1+fso(35):fso(2):ns));
  102.    na=length(sp);
  103.    fso(35)=fso(35)+na*fso(2)-ns;
  104.    if (na)
  105.       fso(31:33)=fso(31:33)+[na;sum(sp);sp.'*sp];
  106.       [sp,fso(36:37)]=filter(1,fso(4:6),sp,fso(36:37));
  107.       % determine new peak
  108.       th=max(sp);
  109.       s=zeros(nb,1);
  110.       h=s;
  111.       if fso(34)>0
  112.          if th>fso(34)
  113.             nt=floor(log(th/fso(34))/log(thf));
  114.             th=fso(34)*thf^nt;
  115.             if nt>-nb
  116.                s(1-nt:nb)=fso(48:47+nb+nt);
  117.                h(1-nt:nb)=fso(56:55+nb+nt);
  118.             end
  119.          else
  120.             th=fso(34);
  121.             s=fso(48:55);
  122.             h=fso(56:63);
  123.          end
  124.       end
  125.       fso(34)=th;
  126.       nh=fso(3);
  127.       if nh<na
  128.          nm=na-nh+2;
  129.          for i=1:nb
  130.             th=th*thf;
  131.             b=sp>th;
  132.             np=h(i);
  133.             s(i)=s(i)+np+sum(cumsum([sum(b(1:np+1));b(2+np:na)]-[zeros(nh-np,1);b(1:na-nh)])>0);
  134.             ff=find([1;b(nm:na)]);
  135.             nm=nm-1+ff(end);
  136.             h(i)=nm+nh-na-2;;
  137.          end
  138.       else
  139.          for i=1:nb
  140.             th=th*thf;
  141.             b=sp>th;
  142.             np=h(i);
  143.             ff=find(b);
  144.             if length(ff)
  145.                h(i)=ff(end)+nh-na-1;
  146.                s(i)=s(i)+na+min(0,np+1-ff(1));
  147.             else
  148.                h(i)=max(0,np-na);
  149.                s(i)=s(i)-h(i)+np;
  150.             end
  151.          end
  152.       end
  153.       fso(48:63)=[s;h];
  154.    end
  155. end
  156. ssq=fso(33);
  157. if ssq>0
  158.    s=fso(48:55);
  159.    thresh=(fso(34)*thf.^(1:nb).').^2;
  160.    % 15.9 dB = 38.9 power ratio
  161.    mar=38.9;
  162.    st=s.*thresh*mar/ssq;
  163.    ff=find(st>1);
  164.    if length(ff)
  165.       nf=ff(end);
  166.       ls=log(s(nf)/s(min(nb,nf+1)));
  167.       lev(1)=ssq/s(nf)*st(nf)^(ls/(ls-2*log(thf)));
  168.    else
  169.       [ls,nf]=max(st);
  170.       lev(1)=ssq/s(nf);
  171.    end
  172.    if length(lev)>1 lev(2)=ssq/fso(31); end
  173. end
  174. if any(mode=='d')
  175.    lev=10*log10(max(lev,1e-200));
  176. end