sigma

   Singularity in EGG by Multiscale Analysis (SIGMA) Algorithm

   [gci goi] = 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, 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 
 Date:             29 May 2007
 Version: $Id: sigma.m 713 2011-10-16 14:45:43Z dmb $
**************************************************************************

0001 function [gci goi] = sigma(lx,fs,fmax)
0002 
0003 %   Singularity in EGG by Multiscale Analysis (SIGMA) Algorithm
0004 %
0005 %   [gci goi] = 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, 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
0035 % Date:             29 May 2007
0036 % Version: $Id: sigma.m 713 2011-10-16 14:45:43Z dmb $
0037 %**************************************************************************
0038 
0039 if(nargin<3)
0040     fmax = 400;
0041 end
0042 
0043 % PARAMETERS - we don't like these
0044 fmin = 50;      % GOI post-filtering
0045 Tmin = 1/fmax;  % Sets GD window length
0046 Tmax = 1/fmin;  % GOI post-filtering
0047 oqmin = 0.1;    % GOI post-filtering
0048 oqmax = 0.9;    % GOI post-filtering
0049 
0050 gwlen = Tmin;   % group delay evaluation window length (normally 0.0025)
0051 fwlen=0.000;    % window length used to smooth group delay (normally 0.000)
0052 
0053 % Normalise (for plotting purposes)
0054 lx = lx/max(abs(lx));
0055 
0056 % Calculate SWT
0057 nlev = 3;
0058 nu = length(lx);
0059 nU=(2^nlev)*ceil(nu./(2^nlev));
0060 [Lo_D Hi_D] = wfilters('bior1.5','d');
0061 [swa swd] = swt([lx; zeros(nU-nu,1)],nlev,Lo_D,Hi_D);
0062 swa=swa(:,1:nu); swd=swd(:,1:nu);
0063 mp = prod(swd)';
0064 mp = [0;mp];
0065 
0066 % Find third roots
0067 nmp = mp;
0068 nmp(find(nmp>0)) = 0;   % Half-wave rectify on negative half of mp for GCI
0069 crnmp = nthroot(nmp,3);
0070 
0071 pmp = mp;
0072 pmp(find(pmp<0)) = 0;   % Half-wave rectify on positive half of mp for GOI
0073 crpmp = nthroot(pmp,3);
0074 
0075 % Group delay evaluation on mp
0076 [gcic sew ngrdel ntoff] = xewgrdel(nmp,fs,gwlen,fwlen);
0077 ngrdel=-[zeros(ntoff,1); ngrdel(1:end-ntoff)];
0078 [goic sew pgrdel ptoff] = xewgrdel(pmp,fs,gwlen,fwlen);
0079 pgrdel=-[zeros(ptoff,1); pgrdel(1:end-ptoff)];
0080 
0081 % Set up other variables
0082 gci = zeros(1,length(lx));      
0083 gciall = zeros(1,length(lx));
0084 gciall(:) = NaN;
0085 gciall(round(gcic)) = 0.25;
0086 
0087 goi = zeros(1,length(lx));      
0088 goiall = zeros(1,length(lx));
0089 goiall(:) = NaN;
0090 goiall(round(goic)) = 0.25;
0091 
0092 % --------- GCI Detection ---------- %
0093 
0094 % Model GD slope
0095 mngrdel = (-((gwlen*fs)/2-1):(gwlen*fs)/2-1)';
0096 mngrdellen = length(mngrdel);
0097 cmngrdel = zeros(length(ngrdel),1);
0098 
0099 for i=1:length(gcic)
0100     lbnd = round(gcic(i)-((gwlen*fs)/2-1)); 
0101     ubnd = lbnd+mngrdellen-1;
0102     
0103     if ~( (lbnd<1) || (ubnd>length(ngrdel)) )
0104         nfv(i,1) = sum(crnmp(lbnd:ubnd));                     % Sum of crnmp over GD window
0105         nfv(i,2) = min(crnmp(lbnd:ubnd));                     % Peak value of crnmp
0106         nfv(i,3) = sqrt( mean( (mngrdel-ngrdel(lbnd:ubnd)).^2 ) ); % Phase slope deviation
0107         
0108         cmngrdel(lbnd:ubnd) = mngrdel;
0109     end
0110 end
0111 
0112 nclust = 2;
0113 snfv = size(nfv);
0114 
0115 % Determine clusters
0116 [mm,vv,ww]=gaussmix(nfv,var(nfv)/100,[],nclust,'kf');
0117 
0118 % Find cluster with lowest crnmp sum
0119 [y I] = min(mm(:,1));
0120 
0121 % Find log likelihoods for each cluster
0122 for i=1:nclust
0123     [m_,v_,w_,g,f,ll(i,:)]=gaussmix(nfv,[],0,mm(i,:),vv(i,:),ww(i,:));
0124 end
0125 [m,in]=max(ll); % Find which cluster each feature vector belongs to
0126 
0127 close all;
0128 naccept = [];
0129 nreject = [];
0130 for i=1:snfv(1)
0131     if (in(i) == I)
0132         naccept = [naccept; nfv(i,:)];
0133     else
0134         nreject = [nreject; nfv(i,:)];
0135     end
0136 end
0137 
0138 % If the candidate belongs to the chosen cluster then keep
0139 for i=1:length(nfv)
0140     if (in(i) == I)
0141         gci(round(gcic(i))) = 0.5;
0142     end
0143 end
0144 
0145 % -------- Post-filter swallows (GCIs only) -------- %
0146 if(length(gci)>2)
0147     % If a gci is separated from all others by more than Tmax, delete.
0148     fgci = find(gci);
0149     % Check first one
0150     if ( (fgci(2)-fgci(1)) > Tmax*fs)
0151         fgci = fgci(2:end);
0152     end
0153     % Check the middle
0154     i=2;
0155     while(i<length(fgci)-1)
0156         if ( ((fgci(i)-fgci(i-1))>Tmax*fs) && ((fgci(i+1)-fgci(i))>Tmax*fs) )
0157             fgci = [fgci(1:i-1) fgci(i+1:end)];
0158         end
0159         i = i+1;
0160     end
0161     % Check last one
0162     if ( (fgci(end)-fgci(end-1)) > Tmax*fs)
0163         fgci = fgci(1:end-1);
0164     end
0165     % Convert back
0166     gci = zeros(1,max(fgci));
0167     gci(fgci) = 0.5;
0168 end
0169 
0170 % --------- GOI Detection ---------- %
0171 % Model GD slope
0172 mpgrdel = (-((gwlen*fs)/2-1):(gwlen*fs)/2-1)';
0173 mpgrdellen = length(mpgrdel);
0174 cmpgrdel = zeros(length(pgrdel),1);
0175 
0176 for i=1:length(goic)
0177     lbnd = round(goic(i)-((gwlen*fs)/2-1)); 
0178     %ubnd = round(goic(i)+((gwlen*fs)/2-1));
0179     ubnd = min(lbnd+mpgrdellen-1,length(crpmp));
0180     
0181     if ~( (lbnd<1) || (ubnd>length(pgrdel)) )
0182         pfv(i,1) = sum(crpmp(lbnd:ubnd));                     % Sum of crnmp over GD window
0183         pfv(i,2) = max(crpmp(lbnd:ubnd));                     % Peak value of crpmp
0184         pfv(i,3) = sqrt( mean( (mpgrdel-pgrdel(lbnd:ubnd)).^2 ) ); % Phase slope deviation
0185                 
0186         cmpgrdel(lbnd:ubnd) = mpgrdel;
0187     end
0188 end
0189 
0190 nclust = 2;
0191 spfv = size(pfv);
0192 
0193 % Determine clusters
0194 [mm,vv,ww]=gaussmix(pfv,var(pfv)/100,[],nclust,'kf');
0195 
0196 % Find cluster with highest crpmp sum
0197 [y I] = max(mm(:,1));
0198 
0199 % Find log likelihoods for each cluster
0200 ll = [];
0201 for i=1:nclust
0202     [m_,v_,w_,g,f,ll(i,:)]=gaussmix(pfv,[],0,mm(i,:),vv(i,:),ww(i,:));
0203 end
0204 [m,in]=max(ll); % Find which cluster each feature vector belongs to
0205 
0206 paccept = [];
0207 preject = [];
0208 for i=1:spfv(1)
0209     if (in(i) == I)
0210         paccept = [paccept; pfv(i,:)];
0211     else
0212         preject = [preject; pfv(i,:)];
0213     end
0214 end
0215 
0216 % If the candidate belongs to the chosen cluster then keep
0217 for i=1:length(pfv)
0218     if (in(i) == I)
0219         goi(round(goic(i))) = 0.75;
0220     end
0221 end
0222 
0223 % ------- GOI Post-Filtering ------- %
0224 
0225 % For all GCIs, find GOIs which are within OQ limits
0226 goiprefilt = goi;
0227 goifilt = zeros(size(goi));
0228 gciind = find(gci);
0229 Tprev = Tmax;
0230 prev = 0;
0231 nofirst = 0;
0232 for i=2:length(gciind)
0233     lbnd = gciind(i-1);
0234     ubnd = gciind(i);
0235     T = ubnd-lbnd;
0236     if(T>Tmax*fs)   % If period is too long, use previous.
0237         T = Tprev;
0238     end
0239     I = find(goi( round(lbnd+(1-oqmax)*T):round(lbnd+(1-oqmin)*T) ));
0240     if(~isempty(I)) % If we have a GOI
0241         prev = round(I(1)+(1-oqmax)*T-1);
0242         goifilt(round(I(1)+lbnd+(1-oqmax)*T-1)) = 0.5;  % Taking first - should it be highest energy?
0243     else            % If not then insert at last OQ.
0244         if(prev>0)
0245             if( (lbnd+prev) < gciind(i) ) % Protect against GOI occuring after next GCI
0246                 goifilt(lbnd + prev-1) = 0.5;
0247             else
0248                 goifilt( gciind(i)-1) = 0.5;
0249             end
0250         end
0251         if(i==2)
0252             nofirst = 1;
0253         end
0254     end
0255     if(nofirst && (prev>0)) % If no GOI occurs after GOI, no previous period exists, so add after a period has been found.
0256         goifilt(gciind(1)+prev-1) = 0.5;
0257         nofirst = 0;
0258     end
0259     Tprev = T;
0260 end
0261 % Final period
0262 lbnd = gciind(end);
0263 I = find(goi( round(lbnd+(1-oqmax)*T):min(round(lbnd+(1-oqmin)*T),length(goi))));
0264 if(~isempty(I))
0265     goifilt(round(I(1)+lbnd+(1-oqmax)*T-1)) = 0.5;
0266 else            % If not then insert at last OQ.
0267     if(prev>0)
0268         goifilt(lbnd + prev) = 0.5;
0269     end
0270 end
0271 goi = goifilt;
0272 
0273 gci = find(gci>0);
0274 goi = find(goi>0);
0275 
0276 %% EW group delay epoch extraction
0277 function [tew,sew,y,toff]=xewgrdel(u,fs,gwlen,fwlen)
0278 
0279 error(nargchk(2,4,nargin));
0280 
0281 if(nargin < 4)
0282     fwlen = voicebox('dy_fwlen');          % window length to smooth group delay
0283 end
0284 if (nargin < 3)
0285     gwlen = voicebox('dy_gwlen');          % window length of group delay
0286 end
0287 
0288 % perform group delay calculation
0289 
0290 gw=2*floor(gwlen*fs/2)+1;            % force window length to be odd
0291 ghw=window('hamming',gw,'s');
0292 ghw = ghw(:);                           % force to be a column (dmb thinks window gives a row - and he should know as he wrote it!)
0293 ghwn=ghw'.*(gw-1:-2:1-gw)/2;            % weighted window: zero in middle
0294 
0295 u2=u.^2;
0296 yn=filter(ghwn,1,u2);
0297 yd=filter(ghw,1,u2);
0298 yd(abs(yd)<eps)=10*eps;                 % prevent infinities
0299 y=yn(gw:end)./yd(gw:end);               % delete filter startup transient
0300 toff=(gw-1)/2;
0301 fw=2*floor(fwlen*fs/2)+1;            % force window length to be odd
0302 if fw>1
0303     daw=window('hamming',fw,'s');
0304     y=filter(daw,1,y)/sum(daw);         % low pass filter
0305     toff=toff-(fw-1)/2;
0306 end
0307 [tew,sew]=zerocros(y,'n');              % find zero crossings
0308 
0309 tew=tew+toff;                           % compensate for filter del