V_LPCRF2AR Convert reflection coefs to autoregressive coefs [AR,ARP,ARU,G]=(RF) Input: RF(:,p+1) gives reflection coefficients of one or more p-section lossless tubes Ouputs: G is the gain of the all-pole AR filter AR/G is the transfer function from U_in to the glottal input wave, U_g. AR(:,1)=1 always. ARP*K is the transfer function from U_in to the pressure just after the glottis where K = rho*c/Alips: rho = air density 1.23 kg/m^3, c=sound speed 340 m/s, Alips = effective area of free space beyond the lips. ARU is the transfer function from U_in to the total volume velocity through the glottis where U_in=z^(p/2)*U_lips is the time-advanced volume velocity at the lips Energy into the vcal tract is equal to K*filter(ARP,1,Ulips).*filter(ARU,1,Ulips) reverse glottal flows divided by 1-r0 where r0 is the glottal reflection coefficient. The scale factor is included to avoid a zero answer when the glottis is closed giving r0=1. The transfer functions have ar(:,1)=art(:,1)=1 They should both be multiplied by z^(p/2)/prod(1+rf) to correct the absolute gain and to compensate for the delay of p/2 samples along the length of the tube. The energy into the vocal tract is given by ars(speech) * are(speech) Ref: D. M. Brookes and H. P. Loke. "Modelling energy flow in the vocal tract with applications to glottal closure and opening detection." In Proc ICASSP'99, pages 213-216, Mar 1999.
0001 function [ar,arp,aru,g]=v_lpcrf2ar(rf) 0002 %V_LPCRF2AR Convert reflection coefs to autoregressive coefs [AR,ARP,ARU,G]=(RF) 0003 % 0004 % Input: RF(:,p+1) gives reflection coefficients of one or more p-section lossless tubes 0005 % Ouputs: G is the gain of the all-pole AR filter 0006 % AR/G is the transfer function from U_in to the glottal input wave, U_g. 0007 % AR(:,1)=1 always. 0008 % ARP*K is the transfer function from U_in to the pressure just after the glottis 0009 % where K = rho*c/Alips: rho = air density 1.23 kg/m^3, c=sound speed 340 m/s, 0010 % Alips = effective area of free space beyond the lips. 0011 % ARU is the transfer function from U_in to the total volume velocity through the glottis 0012 % 0013 % where U_in=z^(p/2)*U_lips is the time-advanced volume velocity at the lips 0014 % 0015 % Energy into the vcal tract is equal to K*filter(ARP,1,Ulips).*filter(ARU,1,Ulips) 0016 % reverse glottal flows divided by 1-r0 where r0 is the glottal reflection coefficient. 0017 % The scale factor is included to avoid a zero answer when the glottis is closed giving r0=1. 0018 % 0019 % The transfer functions have ar(:,1)=art(:,1)=1 0020 % They should both be multiplied by z^(p/2)/prod(1+rf) to correct the absolute 0021 % gain and to compensate for the delay of p/2 samples along the length of the tube. 0022 % 0023 % The energy into the vocal tract is given by ars(speech) * are(speech) 0024 % 0025 % Ref: D. M. Brookes and H. P. Loke. "Modelling energy flow in the vocal tract with 0026 % applications to glottal closure and opening detection." In Proc ICASSP'99, pages 213-216, Mar 1999. 0027 0028 0029 % Copyright (C) Mike Brookes 1997 0030 % Version: $Id: v_lpcrf2ar.m 10865 2018-09-21 17:22:45Z dmb $ 0031 % 0032 % VOICEBOX is a MATLAB toolbox for speech processing. 0033 % Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html 0034 % 0035 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0036 % This program is free software; you can redistribute it and/or modify 0037 % it under the terms of the GNU General Public License as published by 0038 % the Free Software Foundation; either version 2 of the License, or 0039 % (at your option) any later version. 0040 % 0041 % This program is distributed in the hope that it will be useful, 0042 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0043 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0044 % GNU General Public License for more details. 0045 % 0046 % You can obtain a copy of the GNU General Public License from 0047 % http://www.gnu.org/copyleft/gpl.html or by writing to 0048 % Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA. 0049 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0050 0051 [nf,p1]=size(rf); 0052 p2=p1+1; 0053 p=p1-1; 0054 pm=p-1; 0055 arf=[ones(nf,1) zeros(nf,p)]; 0056 arr=[zeros(nf,p) rf(:,p1)]; 0057 cr=zeros(nf,p); 0058 for k=1:p-1 0059 rk=rf(:,(p1-k)*ones(1,k)); 0060 cr(:,1:k)=arr(:,p2-k:p1); 0061 arr(:,p1-k:p)=arr(:,p1-k:p)+rk.*arf(:,1:k); 0062 arf(:,2:k+1)=arf(:,2:k+1)+rk.*cr(:,1:k); 0063 end 0064 r1=rf(:,1); 0065 ar=arf+r1(:,ones(1,p1)).*arr; 0066 if nargout>1 0067 kp=prod(1-rf(:,2:p1),2); 0068 arp=(arf-arr)./kp(:,ones(1,p1)); 0069 if nargout>2 0070 g=prod(1+rf(:,2:p1),2); 0071 aru=(arf+arr)./g(:,ones(1,p1)); 0072 if nargout>3 0073 g=g.*(1+rf(:,1)); 0074 end 0075 end 0076 end 0077 0078