glotlf

GLOTLF   Liljencrants-Fant glottal model U=(D,T,P)
 d is derivative of flow waveform: must be 0, 1 or 2
 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 of 3 parameters defining the waveform
    p(1) is the time at which ugd has its peak negative value. This we define as the
         start of the closed phase. p(1) is therefore the open/closed interval ratio.
    p(2) is the reciprocal of the peak negative value of ugd(t)
    p(3) is the fraction of the open phase for which ugd(t) is negative. That is, it is
         it is the time between the peak flow and the end of the open phase expressed
         as a fraction of the open phase.

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

 Usage example:    ncyc=5;
            period=80;
            t=0:1/period:ncyc;
            ug=glotlf(0,t);
            plot(t,ug)

0001 function u=glotlf(d,t,p)
0002 %GLOTLF   Liljencrants-Fant glottal model U=(D,T,P)
0003 % d is derivative of flow waveform: must be 0, 1 or 2
0004 % t is a vector of time instants at which to calculate the
0005 %   waveform. Time units are in fractions of a cycle.
0006 % p is a vector of 3 parameters defining the waveform
0007 %    p(1) is the time at which ugd has its peak negative value. This we define as the
0008 %         start of the closed phase. p(1) is therefore the open/closed interval ratio.
0009 %    p(2) is the reciprocal of the peak negative value of ugd(t)
0010 %    p(3) is the fraction of the open phase for which ugd(t) is negative. That is, it is
0011 %         it is the time between the peak flow and the end of the open phase expressed
0012 %         as a fraction of the open phase.
0013 %
0014 % Note: this signal has not been low-pass filtered
0015 % and will therefore be aliased [this is a bug]
0016 %
0017 % Usage example:    ncyc=5;
0018 %            period=80;
0019 %            t=0:1/period:ncyc;
0020 %            ug=glotlf(0,t);
0021 %            plot(t,ug)
0022 
0023 
0024 %      Copyright (C) Mike Brookes 1998
0025 %      Version: $Id: glotlf.m 713 2011-10-16 14:45:43Z dmb $
0026 %
0027 %   VOICEBOX is a MATLAB toolbox for speech processing.
0028 %   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
0029 %
0030 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0031 %   This program is free software; you can redistribute it and/or modify
0032 %   it under the terms of the GNU General Public License as published by
0033 %   the Free Software Foundation; either version 2 of the License, or
0034 %   (at your option) any later version.
0035 %
0036 %   This program is distributed in the hope that it will be useful,
0037 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
0038 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0039 %   GNU General Public License for more details.
0040 %
0041 %   You can obtain a copy of the GNU General Public License from
0042 %   http://www.gnu.org/copyleft/gpl.html or by writing to
0043 %   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
0044 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0045 
0046 if nargin < 2
0047   tt=(0:99)'/100;
0048 else
0049   tt=t-floor(t);
0050 end
0051 u=zeros(size(tt));
0052 de=[0.6 0.1 0.2]';
0053 if nargin < 3
0054   p=de;
0055 elseif length(p)<2
0056   p=[p(:); de(length(p)+1:2)];
0057 end
0058 
0059 % Calculate the parameters in terms of ugd(t), the glottal flow derivative
0060 
0061 te=p(1);            % ugd(te) is the negative peak.
0062 mtc=te-1;
0063 e0=1;
0064 wa=pi/(te*(1-p(3)));
0065 a=-log(-p(2)*sin(wa*te))/te;
0066 inta=e0*((wa/tan(wa*te)-a)/p(2)+wa)/(a^2+wa^2);
0067 
0068 % if inta<0 we should reduce p(2)
0069 % if inta>0.5*p(2)*(1-te) we should increase p(2)
0070 
0071 rb0=p(2)*inta;
0072 rb=rb0;
0073 
0074 % Use Newton to determine closure time constant
0075 % so that flow starts and ends at zero.
0076 
0077 for i=1:4
0078   kk=1-exp(mtc/rb);
0079   err=rb+mtc*(1/kk-1)-rb0;
0080   derr=1-(1-kk)*(mtc/rb/kk)^2;
0081   rb=rb-err/derr;
0082 end
0083 e1=1/(p(2)*(1-exp(mtc/rb)));
0084 
0085 
0086 ta=tt<te;
0087 tb=~ta;
0088 
0089 if d==0
0090   u(ta)=e0*(exp(a*tt(ta)).*(a*sin(wa*tt(ta))-wa*cos(wa*tt(ta)))+wa)/(a^2+wa^2);
0091   u(tb)=e1*(exp(mtc/rb)*(tt(tb)-1-rb)+exp((te-tt(tb))/rb)*rb);
0092 elseif d==1
0093   u(ta)=e0*exp(a*tt(ta)).*sin(wa*tt(ta));
0094   u(tb)=e1*(exp(mtc/rb)-exp((te-tt(tb))/rb));
0095 elseif d==2
0096   u(ta)=e0*exp(a*tt(ta)).*(a*sin(wa*tt(ta))+wa*cos(wa*tt(ta)));
0097   u(tb)=e1*exp((te-tt(tb))/rb)/rb;
0098 else
0099   error('Derivative must be 0,1 or 2');
0100 end