activlev.m

function [lev,fso]=activlev(sp,fs,mode)
%ACTIVLEV Measure active speech level as in ITU-T P.56 [LEV,FSO]=(SP,FS,MODE)
%
%Inputs: SP is the speech signal (with better than 20dB SNR)
%        FS is the sample frequency in Hz (see also FSO below)
%        MODE is any combination of the following:
%            r - raw: do not apply the input filter 0.2 to 5.5 kHz
%            a - all: use all the samples rather than subsampling to 694Hz
%            d - give outputs in dB rather than power
%            l - give the "long term level" as well as the active speech level
%Outputs:
%        LEV gives the speech level in units of power (or dB if mode='d')
%            if mode='l' is specified, LEV is a row vector with the "long term
%            level" as its second element (this is just the mean power)
%        FSO is a column vector of intermediate information that allows
%            you to process a speech signal in chunks. Thus:
%
%            fso=fs; for i=1:inc:nsamp, [lev,fso]=activlev(sp(i:i+inc-1),fso,mode); end
%
%            is equivalent to:                lev=activlev(sp(1:nsamp),fs,mode)
%
%            but is much slower. The two methods will not give identical results
%            because they will use slightly different threshods.

%For completeness we list here the contents of the FSO array:
%
%     1 : sample frequency
%     2 : subsampling increment
%     3 : hangover increment
%   4-6 : smoothing filter coefs
%  7-12 : 200Hz HP filter numerator
% 13-18 : 200Hz HP filter denominator
% 19-24 : 5.5kHz LP filter numerator
% 25-30 : 5.5kHz LP filter denominator
% 31-33 : count,sum and sum of squares of filter speech
%    34 : upper limit of largest amplitude bin
%    35 : samples to skip for subsampling
% 36-37 : smoothing filter state
% 38-42 : 200Hz HP filter state
% 43-47 : 5.5kHz LP filter state
% 48-55 : thresholded sample counts
% 56-63 : hangover counts

%      Copyright (C) Mike Brookes 1998
%
%      Last modified Thu Apr 30 17:28:29 1998
%
%   VOICEBOX home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   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<3 mode='0'; else mode = [mode '0']; end
lev=zeros(1,1+any(mode=='l'));
thf=0.5;
nb=8;
fso=[fs;zeros(64-1-length(fs),1)];
if length(fs)==1
   ni=max(floor(fs*all(mode~='a')/694),1);
   ti=ni/fs;
   nh=ceil(0.2/ti)+1;
   g=exp(-ti/0.03);
   % s-plane zeros and poles of high pass 5'th order chebychev2 filter with -0.25dB at w=1
   szp=[0.37843443673309i 0.23388534441447i; -0.20640255179496+0.73942185906851i -0.54036889596392+0.45698784092898i]; 
   szp=[[0; -0.66793268833792] szp conj(szp)];
   zl=2./(1-szp*tan(200*pi/fs))-1;
   al=real(poly(zl(2,:)));
   bl=real(poly(zl(1,:)));
   sw=1-2*rem(0:5,2).';
   bl=bl*(al*sw)/(bl*sw);
   fso(2:19-1)=[ni;nh;[1; -2*g; g^2]/(1-2*g+g^2);bl(:);al(:)];
   if fs>14000
      zh=2./(szp/tan(5500*pi/fs)-1)+1;
      ah=real(poly(zh(2,:)));
      bh=real(poly(zh(1,:)));
      bh=bh*sum(ah)/sum(bh);
      fso(19:19+11)=[bh(:);ah(:)];
   end
end
% need to calculate offset correctly
% subsample the signal
ns=length(sp);
if ns
   % input filter goes here
   if all(mode~='r')
      [sp,fso(38:42)]=filter(fso(7:12),fso(13:18),sp(:),fso(38:42));
      if fso(25)
         [sp,fso(43:47)]=filter(fso(19:24),fso(25:30),sp,fso(43:47));
      end
   end
   sp=abs(sp(1+fso(35):fso(2):ns));
   na=length(sp);
   fso(35)=fso(35)+na*fso(2)-ns;
   if (na)
      fso(31:33)=fso(31:33)+[na;sum(sp);sp.'*sp];
      [sp,fso(36:37)]=filter(1,fso(4:6),sp,fso(36:37));
      % determine new peak
      th=max(sp);
      s=zeros(nb,1);
      h=s;
      if fso(34)>0
         if th>fso(34)
            nt=floor(log(th/fso(34))/log(thf));
            th=fso(34)*thf^nt;
            if nt>-nb
               s(1-nt:nb)=fso(48:47+nb+nt);
               h(1-nt:nb)=fso(56:55+nb+nt);
            end
         else
            th=fso(34);
            s=fso(48:55);
            h=fso(56:63);
         end
      end
      fso(34)=th;
      nh=fso(3);
      if nh<na
         nm=na-nh+2;
         for i=1:nb
            th=th*thf;
            b=sp>th;
            np=h(i);
            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);
            ff=find([1;b(nm:na)]);
            nm=nm-1+ff(end);
            h(i)=nm+nh-na-2;;
         end
      else
         for i=1:nb
            th=th*thf;
            b=sp>th;
            np=h(i);
            ff=find(b);
            if length(ff)
               h(i)=ff(end)+nh-na-1;
               s(i)=s(i)+na+min(0,np+1-ff(1));
            else
               h(i)=max(0,np-na);
               s(i)=s(i)-h(i)+np;
            end
         end
      end
      fso(48:63)=[s;h];
   end
end
ssq=fso(33);
if ssq>0
   s=fso(48:55);
   thresh=(fso(34)*thf.^(1:nb).').^2;
   % 15.9 dB = 38.9 power ratio
   mar=38.9;
   st=s.*thresh*mar/ssq;
   ff=find(st>1);
   if length(ff)
      nf=ff(end);
      ls=log(s(nf)/s(min(nb,nf+1)));
      lev(1)=ssq/s(nf)*st(nf)^(ls/(ls-2*log(thf)));
   else
      [ls,nf]=max(st);
      lev(1)=ssq/s(nf);
   end
   if length(lev)>1 lev(2)=ssq/fso(31); end
end
if any(mode=='d')
   lev=10*log10(max(lev,1e-200));
end