Home > voicebox > v_sigma.m

v_sigma

PURPOSE ^

Singularity in EGG by Multiscale Analysis (SIGMA) Algorithm

SYNOPSIS ^

function [gci goi] = v_sigma(lx,fs,fmax)

DESCRIPTION ^

   Singularity in EGG by Multiscale Analysis (SIGMA) Algorithm

   [gci goi] = v_sigma(lx,fs,fmax)

   Inputs:
       lx      Nx1 vector LX signal
       fs      Sampling freq (Hz)
       fmax    [Optional] max laryngeal freq
   Outputs:
       gci      Vector of gcis as sample instants.
       goi      Vector of gois as sample instants.

   External Functions:
       Function gaussmix in Voicebox.

   Notes:
       Due to the initialization of the EM algorithm on a random data
       point, V_SIGMA does not always produce deterministic behaviour.

   References:
       M. R. P. Thomas and P. A. Naylor, "The SIGMA Algorithm: A Glottal
       Activity Detector for Electroglottographic Signals," IEEE Trans.
       Audio, Speech, Lang. Process., vol.17, no.8, pp.1557-1566, Nov.
       2009.

   Revision History:
       13/09/10: Swallow postfilter threw an error when fewer than 3 GCIs 
       were detected. Final GOI cycle threw an error when search bounds
       exceeded length of signal. Similar problem in line 179 fixed.

**************************************************************************
 Author:           M. R. P. Thomas, 29 May 2007
**************************************************************************
      Copyright (C) M. R. P. Thomas, Mike Brookes 2007-2013
      Version: $Id: v_sigma.m 9702 2017-04-19 10:48:12Z dmb $

   VOICEBOX is a MATLAB toolbox for speech processing.
   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
   http://www.gnu.org/copyleft/gpl.html or by writing to
   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function [gci goi] = v_sigma(lx,fs,fmax)
0002 
0003 %   Singularity in EGG by Multiscale Analysis (SIGMA) Algorithm
0004 %
0005 %   [gci goi] = v_sigma(lx,fs,fmax)
0006 %
0007 %   Inputs:
0008 %       lx      Nx1 vector LX signal
0009 %       fs      Sampling freq (Hz)
0010 %       fmax    [Optional] max laryngeal freq
0011 %   Outputs:
0012 %       gci      Vector of gcis as sample instants.
0013 %       goi      Vector of gois as sample instants.
0014 %
0015 %   External Functions:
0016 %       Function gaussmix in Voicebox.
0017 %
0018 %   Notes:
0019 %       Due to the initialization of the EM algorithm on a random data
0020 %       point, V_SIGMA does not always produce deterministic behaviour.
0021 %
0022 %   References:
0023 %       M. R. P. Thomas and P. A. Naylor, "The SIGMA Algorithm: A Glottal
0024 %       Activity Detector for Electroglottographic Signals," IEEE Trans.
0025 %       Audio, Speech, Lang. Process., vol.17, no.8, pp.1557-1566, Nov.
0026 %       2009.
0027 %
0028 %   Revision History:
0029 %       13/09/10: Swallow postfilter threw an error when fewer than 3 GCIs
0030 %       were detected. Final GOI cycle threw an error when search bounds
0031 %       exceeded length of signal. Similar problem in line 179 fixed.
0032 %
0033 %**************************************************************************
0034 % Author:           M. R. P. Thomas, 29 May 2007
0035 %**************************************************************************
0036 %      Copyright (C) M. R. P. Thomas, Mike Brookes 2007-2013
0037 %      Version: $Id: v_sigma.m 9702 2017-04-19 10:48:12Z dmb $
0038 %
0039 %   VOICEBOX is a MATLAB toolbox for speech processing.
0040 %   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
0041 %
0042 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0043 %   This program is free software; you can redistribute it and/or modify
0044 %   it under the terms of the GNU General Public License as published by
0045 %   the Free Software Foundation; either version 2 of the License, or
0046 %   (at your option) any later version.
0047 %
0048 %   This program is distributed in the hope that it will be useful,
0049 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
0050 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0051 %   GNU General Public License for more details.
0052 %
0053 %   You can obtain a copy of the GNU General Public License from
0054 %   http://www.gnu.org/copyleft/gpl.html or by writing to
0055 %   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
0056 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0057 
0058 if(nargin<3)
0059     fmax = 400;
0060 end
0061 
0062 % PARAMETERS - we don't like these
0063 fmin = 50;      % GOI post-filtering
0064 Tmin = 1/fmax;  % Sets GD window length
0065 Tmax = 1/fmin;  % GOI post-filtering
0066 oqmin = 0.1;    % GOI post-filtering
0067 oqmax = 0.9;    % GOI post-filtering
0068 
0069 gwlen = Tmin;   % group delay evaluation window length (normally 0.0025)
0070 fwlen=0.000;    % window length used to smooth group delay (normally 0.000)
0071 
0072 % Normalise (for plotting purposes)
0073 lx = lx/max(abs(lx));
0074 
0075 % Calculate SWT
0076 nlev = 3;
0077 nu = length(lx);
0078 nU=(2^nlev)*ceil(nu./(2^nlev));
0079 [Lo_D Hi_D] = wfilters('bior1.5','d');
0080 [swa swd] = swt([lx; zeros(nU-nu,1)],nlev,Lo_D,Hi_D);
0081 swa=swa(:,1:nu); swd=swd(:,1:nu);
0082 mp = prod(swd)';
0083 mp = [0;mp];
0084 
0085 % Find third roots
0086 nmp = mp;
0087 nmp(find(nmp>0)) = 0;   % Half-wave rectify on negative half of mp for GCI
0088 crnmp = nthroot(nmp,3);
0089 
0090 pmp = mp;
0091 pmp(find(pmp<0)) = 0;   % Half-wave rectify on positive half of mp for GOI
0092 crpmp = nthroot(pmp,3);
0093 
0094 % Group delay evaluation on mp
0095 [gcic sew ngrdel ntoff] = xewgrdel(nmp,fs,gwlen,fwlen);
0096 ngrdel=-[zeros(ntoff,1); ngrdel(1:end-ntoff)];
0097 [goic sew pgrdel ptoff] = xewgrdel(pmp,fs,gwlen,fwlen);
0098 pgrdel=-[zeros(ptoff,1); pgrdel(1:end-ptoff)];
0099 
0100 % Set up other variables
0101 gci = zeros(1,length(lx));      
0102 gciall = zeros(1,length(lx));
0103 gciall(:) = NaN;
0104 gciall(round(gcic)) = 0.25;
0105 
0106 goi = zeros(1,length(lx));      
0107 goiall = zeros(1,length(lx));
0108 goiall(:) = NaN;
0109 goiall(round(goic)) = 0.25;
0110 
0111 % --------- GCI Detection ---------- %
0112 
0113 % Model GD slope
0114 mngrdel = (-((gwlen*fs)/2-1):(gwlen*fs)/2-1)';
0115 mngrdellen = length(mngrdel);
0116 cmngrdel = zeros(length(ngrdel),1);
0117 
0118 for i=1:length(gcic)
0119     lbnd = round(gcic(i)-((gwlen*fs)/2-1)); 
0120     ubnd = lbnd+mngrdellen-1;
0121     
0122     if ~( (lbnd<1) || (ubnd>length(ngrdel)) )
0123         nfv(i,1) = sum(crnmp(lbnd:ubnd));                     % Sum of crnmp over GD window
0124         nfv(i,2) = min(crnmp(lbnd:ubnd));                     % Peak value of crnmp
0125         nfv(i,3) = sqrt( mean( (mngrdel-ngrdel(lbnd:ubnd)).^2 ) ); % Phase slope deviation
0126         
0127         cmngrdel(lbnd:ubnd) = mngrdel;
0128     end
0129 end
0130 
0131 nclust = 2;
0132 snfv = size(nfv);
0133 
0134 % Determine clusters
0135 [mm,vv,ww]=gaussmix(nfv,var(nfv)/100,[],nclust,'kf');
0136 
0137 % Find cluster with lowest crnmp sum
0138 [y I] = min(mm(:,1));
0139 
0140 % Find log likelihoods for each cluster
0141 for i=1:nclust
0142     [m_,v_,w_,g,f,ll(i,:)]=gaussmix(nfv,[],0,mm(i,:),vv(i,:),ww(i,:));
0143 end
0144 [m,in]=max(ll); % Find which cluster each feature vector belongs to
0145 
0146 % close all; % commented out by dmb, 12/11/2013
0147 naccept = [];
0148 nreject = [];
0149 for i=1:snfv(1)
0150     if (in(i) == I)
0151         naccept = [naccept; nfv(i,:)];
0152     else
0153         nreject = [nreject; nfv(i,:)];
0154     end
0155 end
0156 
0157 % If the candidate belongs to the chosen cluster then keep
0158 for i=1:length(nfv)
0159     if (in(i) == I)
0160         gci(round(gcic(i))) = 0.5;
0161     end
0162 end
0163 
0164 % -------- Post-filter swallows (GCIs only) -------- %
0165 if(length(gci)>2)
0166     % If a gci is separated from all others by more than Tmax, delete.
0167     fgci = find(gci);
0168     % Check first one
0169     if ( (fgci(2)-fgci(1)) > Tmax*fs)
0170         fgci = fgci(2:end);
0171     end
0172     % Check the middle
0173     i=2;
0174     while(i<length(fgci)-1)
0175         if ( ((fgci(i)-fgci(i-1))>Tmax*fs) && ((fgci(i+1)-fgci(i))>Tmax*fs) )
0176             fgci = [fgci(1:i-1) fgci(i+1:end)];
0177         end
0178         i = i+1;
0179     end
0180     % Check last one
0181     if ( (fgci(end)-fgci(end-1)) > Tmax*fs)
0182         fgci = fgci(1:end-1);
0183     end
0184     % Convert back
0185     gci = zeros(1,max(fgci));
0186     gci(fgci) = 0.5;
0187 end
0188 
0189 % --------- GOI Detection ---------- %
0190 % Model GD slope
0191 mpgrdel = (-((gwlen*fs)/2-1):(gwlen*fs)/2-1)';
0192 mpgrdellen = length(mpgrdel);
0193 cmpgrdel = zeros(length(pgrdel),1);
0194 
0195 for i=1:length(goic)
0196     lbnd = round(goic(i)-((gwlen*fs)/2-1)); 
0197     %ubnd = round(goic(i)+((gwlen*fs)/2-1));
0198     ubnd = min(lbnd+mpgrdellen-1,length(crpmp));
0199     
0200     if ~( (lbnd<1) || (ubnd>length(pgrdel)) )
0201         pfv(i,1) = sum(crpmp(lbnd:ubnd));                     % Sum of crnmp over GD window
0202         pfv(i,2) = max(crpmp(lbnd:ubnd));                     % Peak value of crpmp
0203         pfv(i,3) = sqrt( mean( (mpgrdel-pgrdel(lbnd:ubnd)).^2 ) ); % Phase slope deviation
0204                 
0205         cmpgrdel(lbnd:ubnd) = mpgrdel;
0206     end
0207 end
0208 
0209 nclust = 2;
0210 spfv = size(pfv);
0211 
0212 % Determine clusters
0213 [mm,vv,ww]=gaussmix(pfv,var(pfv)/100,[],nclust,'kf');
0214 
0215 % Find cluster with highest crpmp sum
0216 [y I] = max(mm(:,1));
0217 
0218 % Find log likelihoods for each cluster
0219 ll = [];
0220 for i=1:nclust
0221     [m_,v_,w_,g,f,ll(i,:)]=gaussmix(pfv,[],0,mm(i,:),vv(i,:),ww(i,:));
0222 end
0223 [m,in]=max(ll); % Find which cluster each feature vector belongs to
0224 
0225 paccept = [];
0226 preject = [];
0227 for i=1:spfv(1)
0228     if (in(i) == I)
0229         paccept = [paccept; pfv(i,:)];
0230     else
0231         preject = [preject; pfv(i,:)];
0232     end
0233 end
0234 
0235 % If the candidate belongs to the chosen cluster then keep
0236 for i=1:length(pfv)
0237     if (in(i) == I)
0238         goi(round(goic(i))) = 0.75;
0239     end
0240 end
0241 
0242 % ------- GOI Post-Filtering ------- %
0243 
0244 % For all GCIs, find GOIs which are within OQ limits
0245 goiprefilt = goi;
0246 goifilt = zeros(size(goi));
0247 gciind = find(gci);
0248 Tprev = Tmax;
0249 prev = 0;
0250 nofirst = 0;
0251 for i=2:length(gciind)
0252     lbnd = gciind(i-1);
0253     ubnd = gciind(i);
0254     T = ubnd-lbnd;
0255     if(T>Tmax*fs)   % If period is too long, use previous.
0256         T = Tprev;
0257     end
0258     I = find(goi( round(lbnd+(1-oqmax)*T):round(lbnd+(1-oqmin)*T) ));
0259     if(~isempty(I)) % If we have a GOI
0260         prev = round(I(1)+(1-oqmax)*T-1);
0261         goifilt(round(I(1)+lbnd+(1-oqmax)*T-1)) = 0.5;  % Taking first - should it be highest energy?
0262     else            % If not then insert at last OQ.
0263         if(prev>0)
0264             if( (lbnd+prev) < gciind(i) ) % Protect against GOI occuring after next GCI
0265                 goifilt(lbnd + prev-1) = 0.5;
0266             else
0267                 goifilt( gciind(i)-1) = 0.5;
0268             end
0269         end
0270         if(i==2)
0271             nofirst = 1;
0272         end
0273     end
0274     if(nofirst && (prev>0)) % If no GOI occurs after GOI, no previous period exists, so add after a period has been found.
0275         goifilt(gciind(1)+prev-1) = 0.5;
0276         nofirst = 0;
0277     end
0278     Tprev = T;
0279 end
0280 % Final period
0281 lbnd = gciind(end);
0282 I = find(goi( round(lbnd+(1-oqmax)*T):min(round(lbnd+(1-oqmin)*T),length(goi))));
0283 if(~isempty(I))
0284     goifilt(round(I(1)+lbnd+(1-oqmax)*T-1)) = 0.5;
0285 else            % If not then insert at last OQ.
0286     if(prev>0)
0287         goifilt(lbnd + prev) = 0.5;
0288     end
0289 end
0290 goi = goifilt;
0291 
0292 gci = find(gci>0);
0293 goi = find(goi>0);
0294 
0295 %% EW group delay epoch extraction
0296 function [tew,sew,y,toff]=xewgrdel(u,fs,gwlen,fwlen)
0297 
0298 error(nargchk(2,4,nargin));
0299 
0300 if(nargin < 4)
0301     fwlen = voicebox('dy_fwlen');          % window length to smooth group delay
0302 end
0303 if (nargin < 3)
0304     gwlen = voicebox('dy_gwlen');          % window length of group delay
0305 end
0306 
0307 % perform group delay calculation
0308 
0309 gw=2*floor(gwlen*fs/2)+1;            % force window length to be odd
0310 ghw=window('hamming',gw);
0311 ghw = ghw(:);                           % force to be a column (dmb thinks window gives a row - and he should know as he wrote it!)
0312 ghwn=ghw'.*(gw-1:-2:1-gw)/2;            % weighted window: zero in middle
0313 
0314 u2=u.^2;
0315 yn=filter(ghwn,1,u2);
0316 yd=filter(ghw,1,u2);
0317 yd(abs(yd)<eps)=10*eps;                 % prevent infinities
0318 y=yn(gw:end)./yd(gw:end);               % delete filter startup transient
0319 toff=(gw-1)/2;
0320 fw=2*floor(fwlen*fs/2)+1;            % force window length to be odd
0321 if fw>1
0322     daw=window('hamming',fw);
0323     y=filter(daw,1,y)/sum(daw);         % low pass filter
0324     toff=toff-(fw-1)/2;
0325 end
0326 [tew,sew]=zerocros(y,'n');              % find zero crossings
0327 
0328 tew=tew+toff;                           % compensate for filter del

Generated on Tue 19-Sep-2017 12:07:31 by m2html © 2003