GLOTLF Liljencrants-Fant glottal model U=(D,T,P) d is derivative of flow waveform: must be 0, 1 or 2 t is in fractions of a cycle p has one row per output point p(:,1)=open phase [0.6] p(:,2)=+ve/-ve slope ratio [0.1] p(:,3)=closure time constant/closed phase [0.2] Note: this signal has not been low-pass filtered and will therefore be aliased 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 in fractions of a cycle 0005 % p has one row per output point 0006 % p(:,1)=open phase [0.6] 0007 % p(:,2)=+ve/-ve slope ratio [0.1] 0008 % p(:,3)=closure time constant/closed phase [0.2] 0009 % Note: this signal has not been low-pass filtered 0010 % and will therefore be aliased 0011 % 0012 % Usage example: ncyc=5; 0013 % period=80; 0014 % t=0:1/period:ncyc; 0015 % ug=glotlf(0,t); 0016 % plot(t,ug) 0017 0018 0019 % Copyright (C) Mike Brookes 1998 0020 % 0021 % Last modified Thu Apr 30 17:22:00 1998 0022 % 0023 % VOICEBOX home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html 0024 % 0025 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0026 % This program is free software; you can redistribute it and/or modify 0027 % it under the terms of the GNU General Public License as published by 0028 % the Free Software Foundation; either version 2 of the License, or 0029 % (at your option) any later version. 0030 % 0031 % This program is distributed in the hope that it will be useful, 0032 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0033 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0034 % GNU General Public License for more details. 0035 % 0036 % You can obtain a copy of the GNU General Public License from 0037 % ftp://prep.ai.mit.edu/pub/gnu/COPYING-2.0 or by writing to 0038 % Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA. 0039 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0040 0041 if nargin < 2 0042 tt=(0:99)'/100; 0043 else 0044 tt=t-floor(t); 0045 end 0046 u=zeros(size(tt)); 0047 de=[0.6 0.1 0.2]'; 0048 if nargin < 3 0049 p=de; 0050 elseif length(p)<2 0051 p=[p(:); de(length(p)+1:2)]; 0052 end 0053 0054 te=p(1); 0055 mtc=te-1; 0056 e0=1; 0057 wa=pi/(te*(1-p(3))); 0058 a=-log(-p(2)*sin(wa*te))/te; 0059 inta=e0*((wa/tan(wa*te)-a)/p(2)+wa)/(a^2+wa^2); 0060 0061 % if inta<0 we should reduce p(2) 0062 % if inta>0.5*p(2)*(1-te) we should increase p(2) 0063 0064 rb0=p(2)*inta; 0065 rb=rb0; 0066 0067 % Use Newton to determine closure time constant 0068 % so that flow starts and ends at zero. 0069 0070 for i=1:4 0071 kk=1-exp(mtc/rb); 0072 err=rb+mtc*(1/kk-1)-rb0; 0073 derr=1-(1-kk)*(mtc/rb/kk)^2; 0074 rb=rb-err/derr; 0075 end 0076 e1=1/(p(2)*(1-exp(mtc/rb))); 0077 0078 0079 ta=tt<te; 0080 tb=~ta; 0081 0082 if d==0 0083 u(ta)=e0*(exp(a*tt(ta)).*(a*sin(wa*tt(ta))-wa*cos(wa*tt(ta)))+wa)/(a^2+wa^2); 0084 u(tb)=e1*(exp(mtc/rb)*(tt(tb)-1-rb)+exp((te-tt(tb))/rb)*rb); 0085 elseif d==1 0086 u(ta)=e0*exp(a*tt(ta)).*sin(wa*tt(ta)); 0087 u(tb)=e1*(exp(mtc/rb)-exp((te-tt(tb))/rb)); 0088 elseif d==2 0089 u(ta)=e0*exp(a*tt(ta)).*(a*sin(wa*tt(ta))+wa*cos(wa*tt(ta))); 0090 u(tb)=e1*exp((te-tt(tb))/rb)/rb; 0091 else 0092 error('Derivative must be 0,1 or 2'); 0093 end