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 tiem (always 1) [1]
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.
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.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 tiem (always 1) [1] 0032 % 0033 % Values are shown in brackets for the default input parameters: p=[0.6 0.1 0.2] 0034 % Note that the equation for the return phase in Fig. 2 of [1] is wrong; the correct equation is given in (11). 0035 % The value of epsilon in this equation 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 % [1] 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. 0046 % Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html 0047 % 0048 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0049 % This program is free software; you can redistribute it and/or modify 0050 % it under the terms of the GNU General Public License as published by 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 the 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 [1]. ugd(te) is the negative peak. 0083 mtc=te-1; % -1 times duration of return phase 0084 e0=1; % E0 from [1] 0085 wa=pi/(te*(1-p(3))); % omega_g from [1] 0086 a=-log(-p(2)*sin(wa*te))/te; % alpha from [1] 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 [1] 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 [1] 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