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 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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%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 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