# v_glotlf

## PURPOSE V_GLOTLF Liljencrants-Fant glottal model U=(D,T,P)

## SYNOPSIS function [u,q]=v_glotlf(d,t,p)

## DESCRIPTION ```V_GLOTLF   Liljencrants-Fant glottal model U=(D,T,P)
Usage: (1) f0=100;                % frequency in Hz
t=linspace(0,0.1,100); % time axis
u=v_glotlf(0,t*f0)       % glottal waveform using default parameters

(2) v_glotlf(1,[],[0.6 0.2 0.25]); % plot graph if no output argument

Inputs:  d selects which derivative of the flow waveform to calculate; 0<=d<=3 
t is a vector of time instants at which to calculate the
waveform. Time units are in fractions of a cycle.
p is a vector, [te, E0/Ee, 1-tp/te], defining the waveform [0.6 0.1 0.2]
See below for the parameter meanings.
The peak negative value of the flow derivative is Ug'(te) = -Ee
The peak flow occurs at tp.

Outputs:  u the output waveform
q a structure with the following fields taken from Fig 2 of :
q.Up       peak value of Ug occurs at tp            [0.979]
q.E0=1     gain in open phase (always 1)            
q.Ei       peak value of Ug' occurs at ti           [3.57]
q.Ee       min value of Ug' is -Ee and occurs at te 
q.alpha    exponent coefficient in open phase       [4.42]
q.epsilon  exponent coefficient in return phase     [22.4]
q.omega    angular frequency in open phase (=pi/tp) [6.55]
q.t0=0     start time (always 0)                    
q.ti       time of peak in Ug'                      [0.33]
q.tp       time of peak in Ug                       [0.48]
q.te       time of peak negative value in Ug'       [0.6]
q.ta       initial slope of return phase is Ee/ta   [0.045]
q.tc=1     cycle end tiem (always 1)                

Values are shown in brackets for the default input parameters: p=[0.6 0.1 0.2]
Note that for the return phase in Fig. 2 is wrong; the correct equation is (11).
The value of epsilon is chosen to ensure the flow derivative, Ug'(t), integrates to zero.

Bug: this signal has not been low-pass filtered and will therefore be aliased [this is a bug]

References
    G. Fant, J. Liljencrants, and Q. Lin. A four-parameter model of glottal flow. STL-QPSR, 26 (4): 1�13, 1985.```

## CROSS-REFERENCE INFORMATION This function calls:
• v_axisenlarge V_AXISENLARGE - enlarge the axes of a figure (f,h)
• v_texthvc V_TEXTHVC - write text on graph with specified alignment and colour
This function is called by:

## SOURCE CODE ```0001 function [u,q]=v_glotlf(d,t,p)
0002 %V_GLOTLF   Liljencrants-Fant glottal model U=(D,T,P)
0003 %  Usage: (1) f0=100;                % frequency in Hz
0004 %             t=linspace(0,0.1,100); % time axis
0005 %             u=v_glotlf(0,t*f0)       % glottal waveform using default parameters
0006 %
0007 %         (2) v_glotlf(1,[],[0.6 0.2 0.25]); % plot graph if no output argument
0008 %
0009 %  Inputs:  d selects which derivative of the flow waveform to calculate; 0<=d<=3 
0010 %           t is a vector of time instants at which to calculate the
0011 %             waveform. Time units are in fractions of a cycle.
0012 %           p is a vector, [te, E0/Ee, 1-tp/te], defining the waveform [0.6 0.1 0.2]
0013 %             See below for the parameter meanings.
0014 %             The peak negative value of the flow derivative is Ug'(te) = -Ee
0015 %             The peak flow occurs at tp.
0016 %
0017 % Outputs:  u the output waveform
0018 %           q a structure with the following fields taken from Fig 2 of :
0019 %                q.Up       peak value of Ug occurs at tp            [0.979]
0020 %                q.E0=1     gain in open phase (always 1)            
0021 %                q.Ei       peak value of Ug' occurs at ti           [3.57]
0022 %                q.Ee       min value of Ug' is -Ee and occurs at te 
0023 %                q.alpha    exponent coefficient in open phase       [4.42]
0024 %                q.epsilon  exponent coefficient in return phase     [22.4]
0025 %                q.omega    angular frequency in open phase (=pi/tp) [6.55]
0026 %                q.t0=0     start time (always 0)                    
0027 %                q.ti       time of peak in Ug'                      [0.33]
0028 %                q.tp       time of peak in Ug                       [0.48]
0029 %                q.te       time of peak negative value in Ug'       [0.6]
0030 %                q.ta       initial slope of return phase is Ee/ta   [0.045]
0031 %                q.tc=1     cycle end tiem (always 1)                
0032 %
0033 %             Values are shown in brackets for the default input parameters: p=[0.6 0.1 0.2]
0034 %             Note that for the return phase in Fig. 2 is wrong; the correct equation is (11).
0035 %             The value of epsilon is chosen to ensure the flow derivative, Ug'(t), integrates to zero.
0036 %
0037 % Bug: this signal has not been low-pass filtered and will therefore be aliased [this is a bug]
0038 %
0039 % References
0040 %     G. Fant, J. Liljencrants, and Q. Lin. A four-parameter model of glottal flow. STL-QPSR, 26 (4): 1�13, 1985.
0041
0042 %      Copyright (C) Mike Brookes 1998-2017
0043 %      Version: \$Id: v_glotlf.m 10865 2018-09-21 17:22:45Z dmb \$
0044 %
0045 %   VOICEBOX is a MATLAB toolbox for speech processing.
0047 %
0048 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0049 %   This program is free software; you can redistribute it and/or modify
0051 %   the Free Software Foundation; either version 2 of the License, or
0052 %   (at your option) any later version.
0053 %
0054 %   This program is distributed in the hope that it will be useful,
0055 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
0056 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0057 %   GNU General Public License for more details.
0058 %
0059 %   You can obtain a copy of the GNU General Public License from
0060 %   http://www.gnu.org/copyleft/gpl.html or by writing to
0061 %   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
0062 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0063
0064 if nargin<1 || isempty(d)
0065     d=0;                        % default output is Ug(t)
0066 end
0067 if nargin < 2 || isempty(t)
0068     tt=(0:99)'/100;             % default time axis
0069 else
0070     tt=t-floor(t);              % only interested in th fractional part of t
0071 end
0072 u=zeros(size(tt));              % space for output
0073 de=[0.6 0.1 0.2]';              % default parameters
0074 if nargin < 3 || isempty(p)
0075     p=de;
0076 elseif length(p)<2
0077     p=[p(:); de(length(p)+1:2)];
0078 end
0079
0080 % Calculate the parameters in terms of ugd(t), the glottal flow derivative
0081
0082 te=p(1);                        % t_e from . ugd(te) is the negative peak.
0083 mtc=te-1;                       % -1 times duration of return phase
0084 e0=1;                           % E0 from 
0085 wa=pi/(te*(1-p(3)));            % omega_g from 
0086 a=-log(-p(2)*sin(wa*te))/te;    % alpha from 
0087 inta=e0*((wa/tan(wa*te)-a)/p(2)+wa)/(a^2+wa^2); % integral of first portion of waveform (0<t<te)
0088
0089 % if inta<0 we should reduce p(2)
0090 % if inta>0.5*p(2)*(1-te) we should increase p(2)
0091
0092 rb0=p(2)*inta;  % initial time constant neglects the offset; correct if rb<<(1-te)
0093 rb=rb0;         % rb is closure time constant = 1/epsilon in 
0094
0095 % Use Newton to determine closure time constant, rb, that flow starts and ends at zero.
0096
0097 for i=1:4
0098     kk=1-exp(mtc/rb);
0099     err=rb+mtc*(1/kk-1)-rb0;
0100     derr=1-(1-kk)*(mtc/rb/kk)^2;
0101     rb=rb-err/derr;               % rb is closure time constant = 1/epsilon in 
0102 end
0103 e1=1/(p(2)*(1-exp(mtc/rb)));
0104 ta=tt<te;
0105 tb=~ta;
0106 if d==0                                 % output Ug(t)
0107     u(ta)=e0*(exp(a*tt(ta)).*(a*sin(wa*tt(ta))-wa*cos(wa*tt(ta)))+wa)/(a^2+wa^2);
0108     u(tb)=e1*(exp(mtc/rb)*(tt(tb)-1-rb)+exp((te-tt(tb))/rb)*rb);
0109     ylab='u(t)';
0110 elseif d==1                             % output Ug'(t)
0111     u(ta)=e0*exp(a*tt(ta)).*sin(wa*tt(ta));
0112     u(tb)=e1*(exp(mtc/rb)-exp((te-tt(tb))/rb));
0113     ylab='u''(t)';
0114 elseif d==2                             % output Ug''(t)
0115     u(ta)=e0*exp(a*tt(ta)).*(a*sin(wa*tt(ta))+wa*cos(wa*tt(ta)));
0116     u(tb)=e1*exp((te-tt(tb))/rb)/rb;
0117     ylab='u''''(t)';
0118 else
0119     error('Derivative must be 0,1 or 2');
0120 end
0121 if nargout~=1
0122     ti=(pi+atan(-wa/a))/wa;             % time of peak in Ug'
0123     tp=pi/wa;                           % time of peak in Ug
0124     q.Up=e0*wa*(exp(a*tp)+1)/(a^2+wa^2); % peak value of Ug occurs at tp
0125     q.E0=1;                             % gain in open phase (always 1)
0126     q.Ei=e0*exp(a*ti).*sin(wa*ti);      % peak value of Ug' occurs at ti
0127     q.Ee=1/p(2);                        % note ugd(te)=-Ee
0128     q.alpha=a;                          % exponent coefficient in open phase
0129     q.epsilon=1/rb;                     % exponent coefficient in return phase
0130     q.omega=wa;                         % angular frequency in open phase (=pi/tp)
0131     q.t0=0;                             % start of cycle
0132     q.ti=ti;                            % time of peak in Ug'
0133     q.tp=tp;                            % time of peak in Ug
0134     q.te=te;                            % time of peak negative value in Ug'
0135     q.ta=rb/(p(2)*e1);                  % initial slope of return phase is Ee/ta
0136     q.tc=1;                             % end of cycle
0137 end
0138 if ~nargout
0139     plot(tt,u,'-b',[tt(1) tt(end)],[0 0],':k');
0140     v_axisenlarge([-1 -1.05]);
0141     xlim=get(gca,'xlim');
0142     ylim=get(gca,'ylim');
0143     hold on
0144     plot([1;1]*[q.ti q.tp q.te],ylim,':k');
0145     v_texthvc(q.ti,0,' ti','lbk');
0146     v_texthvc(q.tp,0,' tp','lbk');
0147     if d==0
0148         plot(xlim,[1;1]*q.Up,':k');
0149         v_texthvc(0,q.Up,' Up','ltk');
0150         v_texthvc(te,0,' te','lbk');
0151     elseif d==1
0152         plot(xlim,[1;1]*[q.Ei -q.Ee],':k');
0153         plot([te te+q.ta te+q.ta],[-q.Ee 0 ylim(2)],':k');
0154         v_texthvc(te,0,' te','rtk');
0155         v_texthvc(te+q.ta,0,' te+ta','lbk');
0156         v_texthvc(0,-q.Ee,' �Ee','lbk');
0157         v_texthvc(0,q.Ei,' Ei','ltk');
0158     elseif d==2
0159         plot(xlim,[1;1]*q.Ee/q.ta,':k');
0160         v_texthvc(0,q.Ee/q.ta,' Ee/ta','ltk');
0161         v_texthvc(te,0,' te','lbk');
0162     end
0163     hold off
0164     xlabel('Fraction of period');
0165     ylabel(ylab);
0166 end```