Home > voicebox > glotlf.m

glotlf

PURPOSE ^

GLOTLF Liljencrants-Fant glottal model U=(D,T,P)

SYNOPSIS ^

function [u,q]=glotlf(d,t,p)

DESCRIPTION ^

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=glotlf(0,t*f0)       % glottal waveform using default parameters

         (2) 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.
           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 [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 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
 [1]    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: This function is called by:

SOURCE CODE ^

0001 function [u,q]=glotlf(d,t,p)
0002 %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=glotlf(0,t*f0)       % glottal waveform using default parameters
0006 %
0007 %         (2) 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
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 [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 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 % [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: glotlf.m 10204 2017-10-06 07:27:41Z 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 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 [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     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     texthvc(q.ti,0,' ti','lbk');
0146     texthvc(q.tp,0,' tp','lbk');
0147     if d==0
0148         plot(xlim,[1;1]*q.Up,':k');
0149         texthvc(0,q.Up,' Up','ltk');
0150         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         texthvc(te,0,' te','rtk');
0155         texthvc(te+q.ta,0,' te+ta','lbk');
0156         texthvc(0,-q.Ee,' –Ee','lbk');
0157         texthvc(0,q.Ei,' Ei','ltk');
0158     elseif d==2
0159         plot(xlim,[1;1]*q.Ee/q.ta,':k');
0160         texthvc(0,q.Ee/q.ta,' Ee/ta','ltk');
0161         texthvc(te,0,' te','lbk');
0162     end
0163     hold off
0164     xlabel('Fraction of period');
0165     ylabel(ylab);
0166 end

Generated on Tue 10-Oct-2017 08:30:10 by m2html © 2003