# v_dypsa

## PURPOSE

V_DYPSA Derive glottal closure instances from speech [gci,goi] = (s,fs)

## SYNOPSIS

function [gci,goi] = v_dypsa(s,fs)

## DESCRIPTION

    V_DYPSA   Derive glottal closure instances from speech [gci,goi] = (s,fs)
       Note: Needs to be combined with a voiced-voiceless detector to eliminate
       spurious closures in unvoiced and silent regions.

       Inputs:
       s     is the speech signal
       fs    is the sampling frequncy

       Outputs:
       gci   is a vector of glottal closure sample numbers
       gco   is a vector of glottal opening sample numbers derived from
             an assumed constant closed-phase fraction

       References:
       [1]  P. A. Naylor, A. Kounoudes, J. Gudnason, and M. Brookes, "Estimation of Glottal Closure
            Instants in Voiced Speech using the DYPSA Algorithm," IEEE Trans on Speech and Audio
            Processing, vol. 15, pp. 34-43, Jan. 2007.
       [2]  M. Brookes, P. A. Naylor, and J. Gudnason, "A Quantitative Assessment of Group Delay Methods
            for Identifying Glottal Closures in Voiced Speech," IEEE Trans on Speech & Audio Processing,
            vol. 14, no. 2, pp. 456-466, Mar. 2006.
       [3]  A. Kounoudes, P. A. Naylor, and M. Brookes, "The DYPSA algorithm for estimation of glottal
            closure instants in voiced speech," in Proc ICASSP 2002, vol. 1, Orlando, 2002, pp. 349-352.
       [4]  C. Ma, Y. Kamp, and L. F. Willems, "A Frobenius norm approach to glottal closure detection
            from the speech signal," IEEE Trans. Speech Audio Processing, vol. 2, pp. 258-265, Apr. 1994.
       [5]  A. Kounoudes, "Epoch Estimation for Closed-Phase Analysis of Speech," PhD Thesis,
            Imperial College, 2001.

## CROSS-REFERENCE INFORMATION

This function calls:

- v_lpcauto V_LPCAUTO performs autocorrelation LPC analysis [AR,E,K]=(S,P,T)
- v_lpcifilt V_LPCIFILT Apply inverse filter to speech signal U=(S,AR,T,DC,FADE)
- v_maxfilt V_MAXFILT find max of an exponentially weighted sliding window [Y,K,Y0]=(X,F,nn,D,X0)
- v_voicebox V_VOICEBOX set global parameters for Voicebox functions Y=(FIELD,VAL)
- v_zerocros V_ZEROCROS finds the zeros crossings in a signal [T,S]=(Y,M,X)

This function is called by:

## SUBFUNCTIONS

- function z = psp(g,fs)
- function i = zcr(x, p)
- function zz = zcrp(xx)
- function [frob]=frobfun(sp,p,m,offset)
- function goi=simplegci2goi(gci,pr)
- function [tew,sew,y,toff]=xewgrdel(u,fs)
- function Cfn=fnrg(gcic,frob,fs)
- function gci=dpgci(gcic, s, Ch, fnwav, fs)

## SOURCE CODE

    0001 function [gci,goi] = v_dypsa(s,fs)
    0002 %V_DYPSA   Derive glottal closure instances from speech [gci,goi] = (s,fs)
    0003 %   Note: Needs to be combined with a voiced-voiceless detector to eliminate
    0004 %   spurious closures in unvoiced and silent regions.
    0005 %
    0006 %   Inputs:
    0007 %   s     is the speech signal
    0008 %   fs    is the sampling frequncy
    0009 %
    0010 %   Outputs:
    0011 %   gci   is a vector of glottal closure sample numbers
    0012 %   gco   is a vector of glottal opening sample numbers derived from
    0013 %         an assumed constant closed-phase fraction
    0014 %
    0015 %   References:
    0016 %   [1]  P. A. Naylor, A. Kounoudes, J. Gudnason, and M. Brookes, "Estimation of Glottal Closure
    0017 %        Instants in Voiced Speech using the DYPSA Algorithm," IEEE Trans on Speech and Audio
    0018 %        Processing, vol. 15, pp. 34-43, Jan. 2007.
    0019 %   [2]  M. Brookes, P. A. Naylor, and J. Gudnason, "A Quantitative Assessment of Group Delay Methods
    0020 %        for Identifying Glottal Closures in Voiced Speech," IEEE Trans on Speech & Audio Processing,
    0021 %        vol. 14, no. 2, pp. 456-466, Mar. 2006.
    0022 %   [3]  A. Kounoudes, P. A. Naylor, and M. Brookes, "The DYPSA algorithm for estimation of glottal
    0023 %        closure instants in voiced speech," in Proc ICASSP 2002, vol. 1, Orlando, 2002, pp. 349-352.
    0024 %   [4]  C. Ma, Y. Kamp, and L. F. Willems, "A Frobenius norm approach to glottal closure detection
    0025 %        from the speech signal," IEEE Trans. Speech Audio Processing, vol. 2, pp. 258-265, Apr. 1994.
    0026 %   [5]  A. Kounoudes, "Epoch Estimation for Closed-Phase Analysis of Speech," PhD Thesis,
    0027 %        Imperial College, 2001.
    0028
    0029 % Algorithm Parameters
    0030 %       The following parameters are defined in v_voicebox()
    0031 %
    0032 %   dy_cpfrac=0.3;           % presumed closed phase fraction of larynx cycle
    0033 %   dy_cproj=0.2;            % cost of projected candidate
    0034 %   dy_cspurt=-0.45;         % cost of a talkspurt
    0035 %   dy_dopsp=1;              % Use phase slope projection (1) or not (0)?
    0036 %   dy_ewdly=0.0008;         % window delay for energy cost function term [~ energy peak delay from closure] (sec)
    0037 %   dy_ewlen=0.003;          % window length for energy cost function term (sec)
    0038 %   dy_ewtaper=0.001;        % taper length for energy cost function window (sec)
    0039 %   dy_fwlen=0.00045;        % window length used to smooth group delay (sec)
    0040 %   dy_fxmax=500;            % max larynx frequency (Hz)
    0041 %   dy_fxmin=50;             % min larynx frequency (Hz)
    0042 %   dy_fxminf=60;            % min larynx frequency (Hz) [used for Frobenius norm only]
    0043 %   dy_gwlen=0.0030;         % group delay evaluation window length (sec)
    0044 %   dy_lpcdur=0.020;         % lpc analysis frame length (sec)
    0045 %   dy_lpcn=2;               % lpc additional poles
    0046 %   dy_lpcnf=0.001;          % lpc poles per Hz (1/Hz)
    0047 %   dy_lpcstep=0.010;        % lpc analysis step (sec)
    0048 %   dy_nbest=5;              % Number of NBest paths to keep
    0049 %   dy_preemph=50;           % pre-emphasis filter frequency (Hz) (to avoid preemphasis, make this very large)
    0050 %   dy_spitch=0.2;           % scale factor for pitch deviation cost
    0051 %   dy_wener=0.3;            % DP energy weighting
    0052 %   dy_wpitch=0.5;           % DP pitch weighting
    0053 %   dy_wslope=0.1;           % DP group delay slope weighting
    0054 %   dy_wxcorr=0.8;           % DP cross correlation weighting
    0055 %   dy_xwlen=0.01;           % cross-correlation length for waveform similarity (sec)
    0056
    0057 %   Revision History:
    0058 %
    0059 %   3.0 - 29 Jun 2006  - Rewrote DP function to improve speed
    0060 %   2.6 - 29 Jun 2006  - Tidied up algorithm parameters
    0061 %   2.4 - 10 Jun 2006  - Made into a single file aand put into VOICEBOX
    0062 %   2.3 - 18 Mar 2005  - Removed 4kHz filtering of phase-slope function
    0063 %   2.2 - 05 Oct 2004  -  dpgci uses the slopes returned from xewgrdel
    0064 %                      -  gdwav from speech with fs<9000 is not filtered
    0065 %                      -  Various outputs and inputs of functions have been
    0066 %                         removed since now there is no plotting
    0067 %   1.0 - 30 Jan 2001  - Initial version [5]
    0068
    0069 %   Bugs:
    0070 %         1. Allow the projections only to extend to the end of the larynx cycle
    0071 %         2. Compensate for false pitch period cost at the start of a voicespurt
    0072 %         3. Should include energy and phase-slope costs for the first closure of a voicespurt
    0073 %         4. should delete candidates that are too close to the beginning or end of speech for the cost measures
    0074 %            currently this is 200 samples fixed in the main routine but it should adapt to window lengths of
    0075 %            cross-correlation, lpc and energy measures.
    0076 %         5. Should have an integrated voiced/voiceless detector
    0077 %         6. Allow v_dypsa to be called in chunks for a long speech file
    0078 %         7. Do forward & backward passes to allow for gradient descent and/or discriminative training
    0079 %         8. Glottal opening approximation does not really belong in this function
    0080 %         9. The cross correlation window is asymmetric (and overcomplex) for no very good reason
    0081 %        10. Incorporate -0.5 factor into dy_wxcorr and abolish erroneous (nx2-1)/(nx2-2) factor
    0082 %        11. Add in talkspurt cost at the beginning rather than the end of a spurt (more efficient)
    0083 %        12. Remove qmin>2 condition from voicespurt start detection (DYPSA 2 compatibility) in two places
    0084 %        13. Include energy and phase-slope costs at the start of a voicespurt
    0085 %        14. Single-closure voicespurt are only allowed if nbest=1 (should always be forbidden)
    0086 %        15. Penultimate closure candidate is always acceptd
    0087 %        16. Final element of gcic, Cfn and Ch is unused
    0088 %        17. Needs to cope better with irregular voicing (e.g. creaky voice)
    0089 %        18. Should give non-integer GCI positions for improved accuracy
    0090 %        19. Remove constraint that first voicespurt cannot begin until qrmax after the first candidate
    0091
    0092 %      Copyright (C) Tasos Kounoudes, Jon Gudnason, Patrick Naylor and Mike Brookes 2006
    0093 %      Version: $Id: v_dypsa.m 10865 2018-09-21 17:22:45Z dmb $
    0094 %
    0095 %   VOICEBOX is a MATLAB toolbox for speech processing.
    0096 %   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
    0097 %
    0098 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0099 %   This program is free software; you can redistribute it and/or modify
    0100 %   it under the terms of the GNU Lesser General Public License as published by
    0101 %   the Free Software Foundation; either version 3 of the License, or
    0102 %   (at your option) any later version.
    0103 %
    0104 %   This program is distributed in the hope that it will be useful,
    0105 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
    0106 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    0107 %   GNU Lesser General Public License for more details.
    0108 %
    0109 %   You can obtain a copy of the GNU Lesser General Public License from
    0110 %   https://www.gnu.org/licenses/ .
    0111 %
    0112 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0113
    0114 % Extract algorithm constants from VOICEBOX
    0115
    0116
    0117 dy_preemph=v_voicebox('dy_preemph');
    0118 dy_lpcstep=v_voicebox('dy_lpcstep');
    0119 dy_lpcdur=v_voicebox('dy_lpcdur');
    0120 dy_dopsp=v_voicebox('dy_dopsp');              % Use phase slope projection (1) or not (0)?
    0121 dy_ewtaper=v_voicebox('dy_ewtaper');        % Prediction order of FrobNorm method  in seconds
    0122 dy_ewlen=v_voicebox('dy_ewlen');        % windowlength of FrobNorm method  in seconds
    0123 dy_ewdly=v_voicebox('dy_ewdly');        % shift for assymetric speech shape at start of voiced cycle
    0124 dy_cpfrac=v_voicebox('dy_cpfrac');        % presumed ratio of larynx cycle that is closed
    0125 dy_lpcnf=v_voicebox('dy_lpcnf');          % lpc poles per Hz (1/Hz)
    0126 dy_lpcn=v_voicebox('dy_lpcn');            % lpc additional poles
    0127 dy_xwlen=v_voicebox('dy_xwlen');            % cross-correlation length for waveform similarity (sec)
    0128 dy_fxminf=v_voicebox('dy_fxminf');            % minimum pitch for Frobenius norm calculation
    0129
    0130 lpcord=ceil(fs*dy_lpcnf+dy_lpcn);       % lpc poles
    0131
    0132 %PreEmphasise input speech
    0133 s_used=filter([1 -exp(-2*pi*dy_preemph/fs)],1,s);
    0134
    0135 % perform LPC analysis, AC method with Hamming windowing
    0136 [ar, e, k] = v_lpcauto(s_used,lpcord,floor([dy_lpcstep dy_lpcdur]*fs));
    0137
    0138 if any(any(isinf(ar)))    % if the data is bad and gives infinite prediction coefficients we return with a warning
    0139     warning('No GCIs returned');
    0140     gci=[];
    0141     return;
    0142 end;
    0143
    0144 % compute the prediction residual
    0145 r = v_lpcifilt(s_used,ar,k);
    0146
    0147 % compute the group delay function:  EW method from reference [2] above
    0148 [zcr_cand,sew,gdwav,toff]=xewgrdel(r,fs);
    0149 gdwav=-[zeros(toff,1); gdwav(1:end-toff)];
    0150 zcr_cand=[round(zcr_cand), ones(size(zcr_cand))];   %flag zero crossing candidates with ones
    0151
    0152 sew=0.5+sew';  %the phase slope cost of each candidate
    0153
    0154 pro_cand=[];
    0155 if dy_dopsp ~= 0
    0156     pro_cand = psp(gdwav,fs);
    0157     pro_cand = [pro_cand, zeros(length(pro_cand),1)]; %flag projected candidates with zeros
    0158     sew =      [sew zeros(1,size(pro_cand,1))];      %the phase slope cost of a projected candidate is zero
    0159 end;
    0160
    0161 %Sort the zero crossing and projected candidates together and remove any candidates that
    0162 %are to close to the start or end of the speech signal because the cost functions
    0163 %need room either side.
    0164
    0165 [gcic,sin] = sortrows([zcr_cand; pro_cand],1);
    0166 sew=sew(sin);
    0167 saf=max([200,dy_xwlen*fs/2+1,fs/dy_fxminf]);
    0168 sin=find(and(saf<gcic,gcic<length(gdwav)-saf));
    0169 gcic=gcic(sin,:);
    0170 sew=sew(sin);
    0171
    0172 % compute the frobenious norm function used for a cost in the DP
    0173 fnwav=frobfun(s_used,dy_ewtaper*fs,dy_ewlen*fs,dy_ewdly*fs);
    0174
    0175 %Dynamic programming, picks the most likely candidates based on the pitch consistency, energy etc.
    0176 [gci] = dpgci(gcic, s_used(:), sew, fnwav, fs);
    0177
    0178 %Evaluate goi ... determined to be dy_cpfrac percentage of the larynx cylce away from the last gci
    0179 goi=simplegci2goi(gci,dy_cpfrac);
    0180
    0181 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0182
    0183 function z = psp(g,fs)
    0184 %PSP  Calculates the phase slope projections of the group delay function
    0185 %   Z = PSP(G) computes the
    0186
    0187 %   Revision History:
    0188 %       V1.0 July 12th 2002:
    0189 %            Nov. 6th 2002 : if statement added to remove "last" midpoint
    0190
    0191 g = g(:);
    0192
    0193 gdot = [diff(g);0];
    0194 gdotdot = [diff(gdot);0];
    0195
    0196 % find the turning points  as follows: [tp_number, index_of_tp, min(1) or max(-1), g(index_of_tp)]
    0197 turningPoints = zcr(gdot);
    0198 turningPoints = [[1:length(turningPoints)]', turningPoints, sign(gdotdot(turningPoints)), g(turningPoints)];
    0199
    0200 % useful for debug/plotting
    0201 %tplot = zeros(length(g),1);
    0202 %tplot(turningPoints(:,1)) = turningPoints(:,2);
    0203
    0204 % find any maxima which are < 0
    0205 negmaxima = turningPoints(find(turningPoints(:,3) == -1 & turningPoints(:,4) < 0 & turningPoints(:,1)~=1),:);  %Change 01.05.2003 JG: The first row can't be included
    0206
    0207 % find the midpoint between the preceding min and the negative max
    0208 nmi = negmaxima(:,1);
    0209 midPointIndex = turningPoints(nmi-1,2) + round(0.5*(turningPoints(nmi,2) - turningPoints(nmi-1,2)));
    0210 midPointValue = g(midPointIndex);
    0211
    0212 % project a zero crossing with unit slope
    0213 nz = midPointIndex - round(midPointValue);
    0214
    0215 % find any minima which are > 0
    0216 posminima = turningPoints(find(turningPoints(:,3) == 1 & turningPoints(:,4) > 0),:);
    0217
    0218 % find the midpoint between the positive min and the following max
    0219 pmi = posminima(:,1);
    0220
    0221 %Remove last midpoint if it is the last sample
    0222 if ~isempty(pmi), if pmi(end)==size(turningPoints,1), pmi=pmi(1:end-1); end; end;
    0223
    0224 midPointIndex = turningPoints(pmi,2) + round(0.5*(turningPoints(pmi+1,2) - turningPoints(pmi,2)));
    0225 midPointValue = g(midPointIndex);
    0226
    0227 % project a zero crossing with unit slope
    0228 pz = midPointIndex - round(midPointValue);
    0229
    0230 z = sort([nz;pz]);
    0231
    0232 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0233
    0234 function i = zcr(x, p)
    0235 %ZCR  Finds the indices in a vector to  zero crossings
    0236 %   I = ZCR(X) finds the indices of vector X which are closest to zero-crossings.
    0237 %   I = ZCR(X, P) finds indices for positive-going zeros-crossings for P=1 and
    0238 %   negative-going zero-crossings for P=0.
    0239
    0240 x = x(:);
    0241
    0242 if (nargin==2)
    0243     if (p==0)
    0244         z1 = zcrp(x);   % find positive going zero-crossings
    0245     elseif (p==1)
    0246         z1 = zcrp(-x);  % find negative going zero-crossings
    0247     else
    0248         error('ZCR: invalid input parameter 2: must be 0 or 1');
    0249     end
    0250 else
    0251     z1 = [zcrp(x); zcrp(-x)];
    0252 end
    0253
    0254 % find crossings when x==0 exactly
    0255 z0 = find( (x(1:length(x))==0) & ([x(2:length(x));0] ~= 0));
    0256
    0257 % concatenate and sort the two types of zero-crossings
    0258 i = sort([z0; z1]);
    0259
    0260 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0261
    0262 function zz = zcrp(xx)  %only used in zcr
    0263 % find positive-going zero-crossing
    0264 z1 = find(diff(sign(xx)) == -2);
    0265 % find which out of current sample or next sample is closer to zero
    0266 [m, z2] = min([abs(xx(z1)), abs(xx(z1+1))], [], 2);
    0267 zz =  z1 -1 + z2;
    0268
    0269 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0270
    0271 function [frob]=frobfun(sp,p,m,offset)
    0272
    0273 % [frob]=frobfun(sp,p,m)
    0274 %
    0275 % sp is the speech signal assumed to be preemphasised
    0276 % p  is the prediction order  : recomended to be 1 ms in above paper
    0277 % m  is the window length     : recomended to be 1 ms in above paper
    0278 % offset is shift for assymetric speech shape at start of voiced cycle -
    0279 % default 1.5ms.
    0280 %
    0281 % This function implements the frobenius norm based measure C defined in [4] below.
    0282 % It equals the square of the Frobenius norm of the m by p+1 data matrix divided by p+1
    0283 %
    0284 % Reference:
    0285 %   [4]  C. Ma, Y. Kamp, and L. F. Willems, "A Frobenius norm approach to glottal closure detection
    0286 %        from the speech signal," IEEE Trans. Speech Audio Processing, vol. 2, pp. 258"265, Apr. 1994.
    0287
    0288
    0289 %   Revision History:
    0290 %       V1.0 July 12th 2002:
    0291 %            Nov. 6th 2002 : if statement added to remove "last" midpoint
    0292
    0293 %force p m and offset to be integers
    0294 p=round(p);
    0295 m=round(m);
    0296 offset=round(offset);
    0297
    0298 w=(p+1)*ones(1,m+p);
    0299 w(1:p)=1:p;
    0300 w(m+1:p+m)=p:-1:1;
    0301
    0302 w=w./(p+1);
    0303 frob=filter(w,1,sp.^2);
    0304 frob(1:(round((p+m-1)/2) + offset))=[];
    0305
    0306 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0307
    0308 function goi=simplegci2goi(gci,pr)
    0309
    0310 % Estimate glottal opening instants by assuming a fixed closed-phase fraction
    0311
    0312 gci=round(gci);
    0313 maxpitch=max(medfilt1(diff(gci),7));
    0314
    0315 % calculate opening instants
    0316 for kg=1:length(gci)-1
    0317     goi(kg)=gci(kg)+min(pr*(gci(kg+1)-gci(kg)),pr*maxpitch);
    0318 end;
    0319 kg=kg+1;
    0320 goi(kg)=round(gci(kg)+pr*(gci(kg)-gci(kg-1)));  %use the previous pitch period instead
    0321 goi=round(goi);
    0322
    0323 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0324
    0325 function [tew,sew,y,toff]=xewgrdel(u,fs)
    0326
    0327 % implement EW group delay epoch extraction
    0328
    0329 dy_gwlen=v_voicebox('dy_gwlen');          % group delay evaluation window length
    0330 dy_fwlen=v_voicebox('dy_fwlen');          % window length used to smooth group delay
    0331
    0332 % perform group delay calculation
    0333
    0334 gw=2*floor(dy_gwlen*fs/2)+1;            % force window length to be odd
    0335 ghw=window('hamming',gw,'s');
    0336 ghw = ghw(:);                           % force to be a column (dmb thinks window gives a row - and he should know as he wrote it!)
    0337 ghwn=ghw'.*(gw-1:-2:1-gw)/2;            % weighted window: zero in middle
    0338
    0339 u2=u.^2;
    0340 yn=filter(ghwn,1,u2);
    0341 yd=filter(ghw,1,u2);
    0342 yd(abs(yd)<eps)=10*eps;                 % prevent infinities
    0343 y=yn(gw:end)./yd(gw:end);               % delete filter startup transient
    0344 toff=(gw-1)/2;
    0345 fw=2*floor(dy_fwlen*fs/2)+1;            % force window length to be odd
    0346 if fw>1
    0347     daw=window('hamming',fw,'s');
    0348     y=filter(daw,1,y)/sum(daw);         % low pass filter
    0349     toff=toff-(fw-1)/2;
    0350 end
    0351 [tew,sew]=v_zerocros(y,'n');              % find zero crossings
    0352
    0353 tew=tew+toff;                           % compensate for filter delay and frame advance
    0354
    0355
    0356 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0357
    0358 function Cfn=fnrg(gcic,frob,fs)
    0359
    0360 %Frobenious Energy Cost
    0361
    0362 dy_fxminf=v_voicebox('dy_fxminf');
    0363 frob=frob(:)';
    0364 mm=round(fs/dy_fxminf);
    0365 mfrob=v_maxfilt(frob,1,mm);
    0366 mfrob=[mfrob(floor(mm/2)+1:end) max(frob(end-ceil(mm/2):end))*ones(1,floor(mm/2))];
    0367 rfr=frob./mfrob;
    0368 Cfn=0.5-rfr(round(gcic));
    0369
    0370 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0371
    0372 function gci=dpgci(gcic, s, Ch, fnwav, fs)
    0373
    0374 %DPGCI   Choose the best Glottal Closure Instances with Dynamic Programming
    0375 %   gci=dpgci(gcic, s(:), Ch, fnwav, fs) returns vectors of sample indices corresponding
    0376 %   to the instants of glottal closure in the speech signal s at sampling frequency fs Hz.
    0377 %
    0378 %   Inputs:
    0379 %   gcic    is a matrix whos first column are the glottal closure instance candidates and
    0380 %           the second column is 1 if the corresponding gci is derived from a zero crossing
    0381 %           but zero if the gci is from a a projected zero crossing
    0382 %   s       is the speech signal - MUST be a column vector
    0383 %   Ch      the phase slope cost of every candidate
    0384 %   fnwav   is the frobenious norm function of s
    0385 %   fs      is the sampling frequncy
    0386 %
    0387 %   Outputs:
    0388 %   gci     is a vector of glottal closure instances chosen by the DP
    0389
    0390
    0391
    0392 %   Revision History:
    0393 %   Bugs:  Constants are hardwired but defined in a structure like pv (defined in grpdelpv)
    0394 %
    0395
    0396 % get algorithm parameters from v_voicebox()
    0397
    0398 dy_fxmin=v_voicebox('dy_fxmin');        % min larynx frequency (Hz)
    0399 dy_fxmax=v_voicebox('dy_fxmax');        % min larynx frequency (Hz)
    0400 dy_xwlen=v_voicebox('dy_xwlen');        % cross-correlation length for waveform similarity (sec)
    0401 dy_nbest=v_voicebox('dy_nbest');        % Number of NBest paths to keep
    0402 dy_spitch=v_voicebox('dy_spitch');              % scale factor for pitch deviation cost
    0403 wproj=v_voicebox('dy_cproj');           % cost of projected candidate
    0404 dy_cspurt=v_voicebox('dy_cspurt');           % cost of a talkspurt
    0405 dy_wpitch=v_voicebox('dy_wpitch');           % DP pitch weighting
    0406 dy_wener=v_voicebox('dy_wener');           % DP energy weighting
    0407 dy_wslope=v_voicebox('dy_wslope');           % DP group delay slope weighting
    0408 dy_wxcorr=v_voicebox('dy_wxcorr');           % DP cross correlation weighting
    0409
    0410
    0411
    0412 %Constants
    0413 Ncand=length(gcic);
    0414 sv2i=-(2*dy_spitch^2)^(-1);                 % scale factor for pitch deviation cost
    0415 nxc=ceil(dy_xwlen*fs);                      % cross correlation window length in samples
    0416 % === should delete any gci's that are within this of the end.
    0417 % === and for energy window
    0418
    0419 %Limit the search:
    0420 qrmin=ceil(fs/dy_fxmax);
    0421 qrmax=floor(fs/dy_fxmin);
    0422
    0423 %Cost and tracking r = current, q = previous, p = preprevious
    0424 cost=zeros(Ncand, dy_nbest); cost(:,:)=inf;    %Cost matrix, one row for each candidate
    0425 maxcost=zeros(Ncand,1); maxcost(:,:)=inf;   %Maximum cost in each row
    0426 imaxcost=ones(Ncand,1);                     %Index of maximum cost
    0427
    0428 prev = ones(Ncand, dy_nbest);                  %index of previous, q candidates
    0429 ind = ones(Ncand, dy_nbest);                   %index of p in row q (from prev)
    0430 qbest = [zeros(Ncand,1), ones(Ncand,2)]; % the minimum cost in any previous q [cost,q,i]
    0431
    0432 Cfn=fnrg(gcic(:,1),fnwav,fs);  %Frob.Energy Cost
    0433
    0434 %Add start and end state
    0435 % === should probably delete candidates that are too close to either end of the input
    0436 % === why do we ever need the additional one at the tail end ?
    0437 gcic=[[gcic(1,1)-qrmax-2 0];gcic;[gcic(end,1)+qrmax+2 0]];
    0438 Cfn=[0 Cfn 0];
    0439 Ch = [0 Ch 0];
    0440
    0441 % first do parallelized version
    0442
    0443
    0444 % rather complicated window specification is for compatibility with DYPSA 2
    0445 % === +1 below is for compatibility - probably a bug
    0446 wavix=(-floor(nxc/2):floor(nxc/2)+1)';                 % indexes for segments [nx2,1]
    0447 nx2=length(wavix);
    0448 sqnx2=sqrt(nx2);
    0449 g_cr=dy_wener*Cfn+dy_wslope*Ch+wproj*(1-gcic(:,2))';  % fixed costs
    0450
    0451 g_n=gcic(:,1)';                  % gci sample number [1,Ncand+2]
    0452 g_pr=gcic(:,2)';                 % unprojected flag [1,Ncand+2]
    0453 g_sqm=zeros(1,Ncand+1);         % stores: sqrt(nx2) * mean for speech similarity waveform
    0454 g_sd=zeros(1,Ncand+1);         % stores: 1/(Std deviation * sqrt(nx2)) for speech similarity waveform
    0455 f_pq=zeros((Ncand+1)*dy_nbest,1);   % (q-p) period for each node
    0456 f_c=repmat(Inf,(Ncand+1)*dy_nbest,1);    % cumulative cost for each node - initialise to inf
    0457 f_c(1)=0;                       % initial cost of zero for starting node
    0458 % f_costs=zeros(Ncand*dy_nbest,6);   % === debugging only remember costs of candidate
    0459 f_f=ones((Ncand+1)*dy_nbest,1);    % previous node in path
    0460 f_fb=ones((Ncand+1),1);    % points back to best end-of-spurt node
    0461 fbestc=0;                       % cost of best end-of-spurt node
    0462
    0463 qmin=2;
    0464 for r=2:Ncand+1
    0465 %     if r==86
    0466 %         r;
    0467 %     end
    0468     r_n=g_n(r);             % sample number of r = current candidate
    0469     rix=dy_nbest*(r-1)+(1:dy_nbest);    % index range within node variables
    0470
    0471     % determine the range of feasible q candidates
    0472     qmin0=qmin;
    0473     qmin=find(g_n(qmin0-1:r-1)<r_n-qrmax);      % qmin is the nearest candidate that is >qrmax away
    0474     qmin=qmin(end)+qmin0-1;             % convert to absolute index of first viable candidate
    0475     qmax=find(g_n(qmin-1:r-1)<=r_n-qrmin);      % qmax is the nearest candidate that is >=qrmin away
    0476     qmax=qmax(end)+qmin-2;
    0477
    0478
    0479     % calculate waveform similarity cost measure statistics
    0480
    0481     sr=s(r_n+wavix);        % note s MUST be a column vector so sr is also
    0482     wsum=sum(sr);
    0483     g_sqm(r)=wsum/sqnx2;                % mean * sqrt(nx2)
    0484     g_sd(r)=1/sqrt(sr.'*sr-wsum^2/nx2);   % 1/(Std deviation * sqrt(nx2))
    0485
    0486     % now process the candidates
    0487
    0488     if qmin<=qmax
    0489         qix=qmin:qmax;      % q index
    0490         nq=length(qix);
    0491         % === should integrate the -0.5 into dy_wxcorr
    0492         % === the factor (nx2-1)/(nx2-2) is to compensate for a bug in swsc()
    0493         q_cas=-0.5*(nx2-1)/(nx2-2)*dy_wxcorr*(sum(s(repmat(g_n(qix),nx2,1)+repmat(wavix,1,nq)).*repmat(sr,1,nq),1)-g_sqm(qix)*g_sqm(r)).*g_sd(qix)*g_sd(r);
    0494         % compare: i=35; Ca=swsc(g_n(qix(i)),g_n(r),s,fs); [i qix(i) r  g_n(qix(i)) g_n(r) dy_wxcorr*Ca q_cas(i)]
    0495
    0496         % now calculate pitch deviation cost
    0497
    0498         fix = 1+(qmin-1)*dy_nbest:qmax*dy_nbest;    % node index range
    0499         f_qr=repmat(r_n-g_n(qix),dy_nbest,1);    % (r-p) period for each node
    0500         f_pr=f_qr(:)+f_pq(fix);
    0501         % === could absorb the 2 into sv2i
    0502         f_nx=2-2*f_pr./(f_pr+abs(f_qr(:)-f_pq(fix)));
    0503         f_cp=dy_wpitch*(0.5-exp(sv2i*f_nx.^2));
    0504         % === fudge to match dypsa2.4 - could more efficiently be added
    0505         % === onto the cost of a talkspurt end
    0506         % === should be a v_voicebox parameter anyway
    0507         f_cp(f_pq(fix)==0)=dy_cspurt*dy_wpitch;
    0508
    0509         % now find the N-best paths
    0510
    0511         [r_cnb,nbix]=sort(f_c(fix)+f_cp+reshape(repmat(q_cas,dy_nbest,1),nq*dy_nbest,1));
    0512         f_c(rix)=r_cnb(1:dy_nbest)+g_cr(r);     % costs
    0513         f_f(rix)=nbix(1:dy_nbest)+(qmin-1)*dy_nbest;       % traceback nodes
    0514         f_pq(rix)=f_qr(nbix(1:dy_nbest));       % previous period
    0515         % === f_costs is only for debugging
    0516 %         r;
    0517 %         f_costs(rix,1)=f_c(fix(nbix(1:dy_nbest)));
    0518 %         f_costs(rix,2)=wproj*(1-gcic(r,2));
    0519 %         f_costs(rix,3)=f_cp(nbix(1:dy_nbest));
    0520 %         f_costs(rix,4)=dy_wener*Cfn(r);
    0521 %         f_costs(rix,5)=dy_wslope*Ch(r);
    0522 %         f_costs(rix,6)=reshape(q_cas(1+floor((nbix(1:dy_nbest)-1)/dy_nbest)),dy_nbest,1);
    0523
    0524         % check cost of using this candidate as the start of a new spurt
    0525         % ==== the qmin>2 condition is for compatibility with v_dypsa 2 and
    0526         % prevents any spurts starting until at least qrmax past the first
    0527         % gci. This is probably a bug (see again below)
    0528         iNb=rix(end);
    0529         if (qmin>2) && (f_c(f_fb(qmin-1))+wproj*(1-gcic(r,2))<f_c(iNb))        % compare with worst of Nbest paths
    0530             f_f(iNb)=f_fb(qmin-1);
    0531             % === for now we exclude the energy and phase-slope costs for compatibility with dypsa2
    0532             % === this is probably a bug
    0533             f_c(iNb)=f_c(f_fb(qmin-1))+wproj*(1-gcic(r,2));     % replace worst of the costs with start voicespurt cost
    0534             f_pq(iNb)=0;                    % false pq period
    0535         end
    0536         if f_c(rix(1))<fbestc
    0537             f_fb(r)=rix(1);                          % points to the node with lowest end-of-spurt cost
    0538             % === should compensate for the pitch period cost incurred at the start of the next spurt
    0539             % === note that a node can never be a one-node voicespurt on its own unless dy_nbest=1
    0540             % since the start voices[purt option replaced the worst Nbest cost. This is probably good but
    0541             % is a bit inconsistent.
    0542             fbestc=f_c(rix(1));
    0543         else
    0544             f_fb(r)=f_fb(r-1);
    0545         end
    0546     else            % no viable candidates - must be the start of a voicespurt if anything
    0547         % === for now we exclude the energy and phase-slope costs for compatibility with dypsa2
    0548         % === this is probably a bug
    0549         % ==== the qmin>2 condition is for compatibility with v_dypsa 2 and
    0550         % prevents any spurts starting until at least qrmax past the first
    0551         % gci. This is probably a bug (see again above)
    0552         if (qmin>2)
    0553             f_c(rix(1))=f_c(f_fb(qmin-1))+wproj*(1-gcic(r,2));  % cost of new voicespurt
    0554             f_f(rix)=f_fb(qmin-1);                              % traceback to previous talkspurt end
    0555             f_pq(rix)=0;                                        % previous period
    0556         end
    0557         f_fb(r)=f_fb(r-1);                                  % cannot be the end of a voicespurt
    0558     end
    0559 end
    0560
    0561 % now do the traceback
    0562
    0563 gci = zeros(1,Ncand+1);
    0564
    0565 % === for compatibility with dypsa2, we force the penultimate candidate to be accepted
    0566 % === should be: i=f_fb(Ncand+1) but instead we pick the best of the penultimate candidate
    0567 i=rix(1)-dy_nbest;
    0568 if f_c(i-dy_nbest+1)<f_c(i)     % check if start of a talkspurt
    0569     i=i-dy_nbest+1;
    0570 end
    0571 k=1;
    0572 while i>1
    0573     j=1+floor((i-1)/dy_nbest);          % convert node number to candidate number
    0574     gci(k)=g_n(j);
    0575     i=f_f(i);
    0576     k=k+1;
    0577 end
    0578 gci=gci(k-1:-1:1);           % put into ascending order
    0579
    0580
    0581 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0582
    0583