# 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 [0]
               t is a vector of time instants at which to calculate the waveform. Time units are in
               fractions of a cycle measured from the start of the open phase (the integer part is ignored).
               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 when the flow derivative equals zero.

     Outputs:  u the output waveform
               q a structure with the following fields taken from Fig 2 of reference [1]:
                    q.Up       peak value of Ug occurs at tp            [0.979]
                    q.E0=1     gain in open phase (always 1)            [1]
                    q.Ei       peak value of Ug' occurs at ti           [3.57]
                    q.Ee       min value of Ug' is -Ee and occurs at te [10]
                    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)                    [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 time (always 1)                [1]
                    q.Utc      Integral of Ug' over one cycle           [approx 0]

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

       In the source-tract model of speech production [2, Ch.3], the glottal flow waveform passes through the vocal tract filter
       (typically an all-pole LPC filter) and then a lip-radiation filter (typically a differentiator: R(z)=1-1/z) which converts
       volume flow at the lips to pressure. If the filters are quasi-constant, their order may be interchanged and the speech
       pressure waveform obtained by applying the vocal tract filter to the first derivative of the glottal flow waveform.

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

     References
     [1]    G. Fant, J. Liljencrants, and Q. Lin. A four-parameter model of glottal flow. STL-QPSR, 26 (4): 1-13, 1985.
     [2]   L. R. Rabiner and R. W. Schafer. Digital Processing of Speech Signals. Prentice-Hall, 1978. ISBN 0-13-213603-1.

## 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 [0]
    0010 %           t is a vector of time instants at which to calculate the waveform. Time units are in
    0011 %           fractions of a cycle measured from the start of the open phase (the integer part is ignored).
    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 when the flow derivative equals zero.
    0016 %
    0017 % Outputs:  u the output waveform
    0018 %           q a structure with the following fields taken from Fig 2 of reference [1]:
    0019 %                q.Up       peak value of Ug occurs at tp            [0.979]
    0020 %                q.E0=1     gain in open phase (always 1)            [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 [10]
    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)                    [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 time (always 1)                [1]
    0032 %                q.Utc      Integral of Ug' over one cycle           [approx 0]
    0033 %
    0034 %             Values are shown in brackets for the default input parameters: p=[0.6 0.1 0.2]
    0035 %             Note that the equation for the return phase in Fig. 2 of [1] is wrong; the correct equation is given in (11).
    0036 %             The value of epsilon in this equation is chosen to ensure the flow derivative, Ug'(t), integrates to zero.
    0037 %
    0038 %   In the source-tract model of speech production [2, Ch.3], the glottal flow waveform passes through the vocal tract filter
    0039 %   (typically an all-pole LPC filter) and then a lip-radiation filter (typically a differentiator: R(z)=1-1/z) which converts
    0040 %   volume flow at the lips to pressure. If the filters are quasi-constant, their order may be interchanged and the speech
    0041 %   pressure waveform obtained by applying the vocal tract filter to the first derivative of the glottal flow waveform.
    0042 %
    0043 % Bug: this signal has not been low-pass filtered and will therefore be aliased [this is a bug]
    0044 %
    0045 % References
    0046 % [1]    G. Fant, J. Liljencrants, and Q. Lin. A four-parameter model of glottal flow. STL-QPSR, 26 (4): 1-13, 1985.
    0047 % [2]   L. R. Rabiner and R. W. Schafer. Digital Processing of Speech Signals. Prentice-Hall, 1978. ISBN 0-13-213603-1.
    0048
    0049 %      Copyright (C) Mike Brookes 1998-2017
    0050 %      Version: $Id: v_glotlf.m 10865 2018-09-21 17:22:45Z dmb $
    0051 %
    0052 %   VOICEBOX is a MATLAB toolbox for speech processing.
    0053 %   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
    0054 %
    0055 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0056 %   This program is free software; you can redistribute it and/or modify
    0057 %   it under the terms of the GNU Lesser General Public License as published by
    0058 %   the Free Software Foundation; either version 3 of the License, or
    0059 %   (at your option) any later version.
    0060 %
    0061 %   This program is distributed in the hope that it will be useful,
    0062 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
    0063 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    0064 %   GNU Lesser General Public License for more details.
    0065 %
    0066 %   You can obtain a copy of the GNU Lesser General Public License from
    0067 %   https://www.gnu.org/licenses/ .
    0068 %
    0069 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0070
    0071 if nargin<1 || isempty(d)
    0072     d=0;                        % default output is Ug(t)
    0073 end
    0074 if nargin < 2 || isempty(t)
    0075     tt=(0:99)'/100;             % default time axis
    0076 else
    0077     tt=t-floor(t);              % only interested in the fractional part of t
    0078 end
    0079 u=zeros(size(tt));              % space for output
    0080 de=[0.6 0.1 0.2]';              % default parameters
    0081 if nargin < 3 || isempty(p)
    0082     p=de;
    0083 elseif length(p)<2
    0084     p=[p(:); de(length(p)+1:2)];
    0085 end
    0086
    0087 % Calculate the parameters in terms of ugd(t), the glottal flow derivative
    0088
    0089 te=p(1);                                % t_e from [1]. ugd(te) is the negative peak.
    0090 mtc=te-1;                               % -1 times duration of return phase
    0091 e0=1;                                   % E0 from [1]
    0092 wa=pi/(te*(1-p(3)));                    % omega_g from [1]
    0093 a=-log(-p(2)*sin(wa*te))/te;            % alpha from [1]
    0094 inta=e0*((wa/tan(wa*te)-a)/p(2)+wa)/(a^2+wa^2); % integral of first portion of waveform (0<t<te)
    0095
    0096 % convergence is only possible if 0 <= inta <= 0.5*(1-te)/p(2) as ta varies between 0 and 1-te
    0097 % if inta<0 we must reduce p(2); if inta>0.5*(1-te)/p(2) we must increase p(2)
    0098
    0099 rb0=p(2)*inta;                          % initial time constant neglects the offset; correct if rb<<(1-te)
    0100 rb=rb0;                                 % rb is closure time constant = 1/epsilon in [1]
    0101
    0102 % Use Newton to determine closure time constant, rb, so that flow starts and ends at zero.
    0103 thresh=1e-9;                            % convergence threshold
    0104 for i=1:15                              % maximum of 15 iterations (usually fewer)
    0105     kk=1-exp(mtc/rb);                   % kk = ta/rb = ta * epsilon in [1]
    0106     err=rb+mtc*(1/kk-1)-rb0;
    0107     derr=1-(1-kk)*(mtc/rb/kk)^2;
    0108     rb=rb-err/derr;                     % rb is closure time constant = 1/epsilon in [1]
    0109     if abs(err)<thresh, break, end
    0110 end
    0111 if abs(err)>thresh                      % print error if unable to find a value of rb that gives a zero integral for U'(t)
    0112     error('Requested glottal waveform parameters are not feasible');
    0113 end
    0114 e1=1/(p(2)*(1-exp(mtc/rb)));
    0115 ta=tt<te;
    0116 tb=~ta;
    0117 if d==0                                 % output Ug(t)
    0118     u(ta)=e0*(exp(a*tt(ta)).*(a*sin(wa*tt(ta))-wa*cos(wa*tt(ta)))+wa)/(a^2+wa^2);
    0119     u(tb)=e1*(exp(mtc/rb)*(tt(tb)-1-rb)+exp((te-tt(tb))/rb)*rb);
    0120     ylab='u(t)';
    0121 elseif d==1                             % output Ug'(t)
    0122     u(ta)=e0*exp(a*tt(ta)).*sin(wa*tt(ta));
    0123     u(tb)=e1*(exp(mtc/rb)-exp((te-tt(tb))/rb));
    0124     ylab='u''(t)';
    0125 elseif d==2                             % output Ug''(t)
    0126     u(ta)=e0*exp(a*tt(ta)).*(a*sin(wa*tt(ta))+wa*cos(wa*tt(ta)));
    0127     u(tb)=e1*exp((te-tt(tb))/rb)/rb;
    0128     ylab='u''''(t)';
    0129 else
    0130     error('Derivative must be 0,1 or 2');
    0131 end
    0132 if nargout~=1
    0133     ti=(pi+atan(-wa/a))/wa;             % time of peak in Ug'
    0134     tp=pi/wa;                           % time of peak in Ug
    0135     q.Up=e0*wa*(exp(a*tp)+1)/(a^2+wa^2); % peak value of Ug occurs at tp
    0136     q.E0=1;                             % gain in open phase (always 1)
    0137     q.Ei=e0*exp(a*ti).*sin(wa*ti);      % max value of Ug': Ug'(ti)=Ei
    0138     q.Ee=1/p(2);                        % negated min value of Ug': Ug'(te)=-Ee
    0139     q.alpha=a;                          % exponent coefficient in open phase
    0140     q.epsilon=1/rb;                     % exponent coefficient in return phase
    0141     q.omega=wa;                         % angular frequency in open phase (omega = pi/tp)
    0142     q.t0=0;                             % start of cycle
    0143     q.ti=ti;                            % time of peak in Ug': Ug'(ti)=Ei
    0144     q.tp=tp;                            % time of peak in Ug : Ug(tp)=Up
    0145     q.te=te;                            % time of peak negative value in Ug': Ug'(te)=-Ee
    0146     q.ta=rb/(p(2)*e1);                  % initial slope of return phase is Ug''(te)=-Ee/ta
    0147     q.tc=1;                             % end of cycle
    0148     q.Utc=-err/p(2);                    % integral of Ug' (should always be zero)
    0149 end
    0150 if ~nargout
    0151     plot(t,u,'-b',[t(1) t(end)],[0 0],':k');
    0152     v_axisenlarge([-1 -1.05]);
    0153     xlim=get(gca,'xlim');
    0154     ylim=get(gca,'ylim');
    0155     hold on
    0156     plot([1;1]*[q.ti q.tp q.te],ylim,':k');
    0157     v_texthvc(q.ti,0,' ti','lbk');
    0158     v_texthvc(q.tp,0,' tp','lbk');
    0159     if d==0
    0160         plot(xlim,[1;1]*q.Up,':k');
    0161         v_texthvc(0,q.Up,' Up','ltk');
    0162         v_texthvc(te,0,' te','lbk');
    0163     elseif d==1
    0164         plot(xlim,[1;1]*[q.Ei -q.Ee],':k');
    0165         plot([te te+q.ta te+q.ta],[-q.Ee 0 ylim(2)],':k');
    0166         v_texthvc(te,0,' te','rtk');
    0167         v_texthvc(te+q.ta,0,' te+ta','lbk');
    0168         v_texthvc(0,-q.Ee,' Ee','lbk');
    0169         v_texthvc(0,q.Ei,' Ei','ltk');
    0170     elseif d==2
    0171         plot(xlim,[1;1]*q.Ee/q.ta,':k');
    0172         v_texthvc(0,q.Ee/q.ta,' Ee/ta','ltk');
    0173         v_texthvc(te,0,' te','lbk');
    0174     end
    0175     hold off
    0176     xlabel('Fraction of period');
    0177     ylabel(ylab);
    0178 end