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.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