# v_sigma

## PURPOSE

V_SIGMA Estimate glottal opening an closing instants

## SYNOPSIS

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

## DESCRIPTION

    V_SIGMA Estimate glottal opening an closing instants
       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 Lesser General Public License as published by
       the Free Software Foundation; either version 3 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 Lesser General Public License for more details.

       You can obtain a copy of the GNU Lesser General Public License from
       https://www.gnu.org/licenses/ .

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

## CROSS-REFERENCE INFORMATION

This function calls:

This function is called by:

## SUBFUNCTIONS

- function [tew,sew,y,toff]=xewgrdel(u,fs,gwlen,fwlen)

## SOURCE CODE

    0001 function [gci goi] = v_sigma(lx,fs,fmax)
    0002 %V_SIGMA Estimate glottal opening an closing instants
    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 Lesser General Public License as published by
    0045 %   the Free Software Foundation; either version 3 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 Lesser General Public License for more details.
    0052 %
    0053 %   You can obtain a copy of the GNU Lesser General Public License from
    0054 %   https://www.gnu.org/licenses/ .
    0055 %
    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