tf.m

这是一个用于语音信号处理的工具箱

function [TF]=tf(area,leng)
% [TF]=tf(area,leng)
% Vocal tract area(1,60) and leng(1,60) are the inputs;
% Output is the complex TF value;

	glotoc = 1;
	
	NING = 60;
	NTS = 29;
	Nfft = 64;
	NTlen = 11;
	NTsno = 11;
	NTslen = 1;
	NTArea(1) = 0.0001;
	NTArea(2) = 2.0;
	NTArea(3) = 3.0;
	NTArea(4) = 4.0;
	NTArea(5) = 6.0;
	NTArea(6) = 8.0;
	NTArea(7) = 8.0;
	NTArea(8) = 6.5;
	NTArea(9) = 4.0;
	NTArea(10) = 2.0;
	NTArea(11)= 2.0;

	NSIN = 0;
	Nsin(1)=0;
	Nsin(2)=0;

	
% function area2RLC(ar, sl)


	w = 2*pi * 4000.0;
	rho = 1.14e-3;
	mu = 1.86e-4;
	c = 3.53e4;
	eta = 1.4;
	lambda = 5.5e-5;
	epsilon = 0.24;
	mass = 1.5;
	viscous = 1500;
	compilance = 0;
	rdem = sqrt(rho*mu*w/2.0);
	rc2 = rho*c*c;
	Gpc = (eta-1.0)*sqrt(lambda*w/(2.0*epsilon*rho))/rc2;

	for j=1:NING,
 
		peri = 4.0 * sqrt(pi*area(j));
		Rs(j) = peri * rdem/((area(j))^2);
		Ls(j) = rho / area(j);
		Cp(j) = area(j) / rc2;
		Gp(j) = peri * Gpc;

		wall = sqrt(pi)*area(j)*leng(j);	% Maeda %
		Rw(j) = viscous / wall;
		Lw(j) = mass / wall;
	end
	Rr = (128.0*rho*c)/(9.0*pi*pi*area(NING));
	Lr = 1.5 * ( 8.0*rho / (3.0*pi*sqrt(pi*area(NING))) );


% function NTarea2RLC()

	w = 2*pi * 4000.0;
	rho = 1.14e-3;
	mu = 1.86e-4;
	c = 3.53e4;
	eta = 1.4;
	lambda = 5.5e-5;
	epsilon = 0.24;
	mass = 1.5;
	viscous = 1500;
	compilance = 0;
	rdem = sqrt(rho*mu*w/2.0);
	rc2 = rho*c*c;
	Gpc = (eta-1.0)*sqrt(lambda*w/(2.0*epsilon*rho))/rc2;

	for j=1:NTsno,
 
		peri = 4.0 * sqrt(pi*NTArea(j));

		RNs(j) = peri * rdem/(NTArea(j)*NTArea(j));
		LNs(j) = rho / NTArea(j);
		CNp(j) = NTArea(j) / rc2;
		GNp(j) = peri * Gpc;

		wall = sqrt(pi)*NTArea(j)*NTslen; % Maeda %
		RNw(j) = viscous / wall;
		LNw(j) = mass / wall;
	end
	RNr = (128.0*rho*c)/(9.0*pi*pi*NTArea(NTsno));
	LNr = 1.5 * ( 8.0*rho / (3.0*pi*sqrt(pi*NTArea(NTsno))) );


% calculate TRANSFER FUNCTION.

df = 2*pi*5000.0/Nfft;
w = df;
l = 1;

while l <= Nfft,

% from velopharyngeal port to lips
		Aoold = 1 + 0*i;
		Boold = 0 + 0*i;
		Coold = 0 + 0*i;
		Doold = 1 + 0*i;
		for j=NTS+1:NING,
			z = Rs(j) + w*Ls(j)*i;
			wa = Rw(j) + w*Lw(j)*i;
			Zwall = Gp(j) + w*Cp(j)*i;
			y = Zwall + 1/wa;
			Z0 = (z/y)^0.5;
			r = (z*y)^0.5;
			Za = Z0*tanh(leng(j)/2.0*r);
			Zb = Z0/sinh(leng(j)*r);
			a11 = 1.0 + Za/Zb;
			a12 = -(2.0*Za + Za*(Za/Zb));
			a21 = -1/Zb;
			a22 = a11;
			Ao = a11*Aoold + a12*Coold;
			Bo = a11*Boold + a12*Doold;
			Co = a21*Aoold + a22*Coold;
			Do = a21*Boold + a22*Doold;
			Aoold = Ao;
			Boold = Bo;
			Coold = Co;
			Doold = Do;
		end
		Zor = (0.0 + w*Rr*Lr*i)/(Rr + w*Lr*i);
		Hnum = Bo - Do*Zor;
		Hdem = Co*Zor - Ao;
		Zo = Hnum / Hdem;

% from velopharyngeal port to nostrils %
		Anold = 1 + 0*i;
		Bnold = 0 + 0*i;
		Cnold = 0 + 0*i;
		Dnold = 1 + 0*i;
		for j=1:NTsno,
			z = RNs(j) + w*LNs(j)*i;
			wa = RNw(j) + w*LNw(j)*i;
			Zwall = GNp(j) + w*CNp(j)*i;
			y = Zwall + 1/wa;
			Z0 = (z/y)^0.5;
			r = (z*y)^0.5;
			Za = Z0*tanh(NTslen/2.0*r);
			Zb = Z0/sinh(NTslen*r);
			a11 = 1.0 + Za/Zb;
			a12 = -(2.0*Za + Za*(Za/Zb));
			a21 = -(1/Zb);
			a22 = a11;
			An = a11*Anold + a12*Cnold;
			Bn = a11*Bnold + a12*Dnold;
			Cn = a21*Anold + a22*Cnold;
			Dn = a21*Bnold + a22*Dnold;
			Anold = An;
			Bnold = Bn;
			Cnold = Cn;
			Dnold = Dn;

% check if any sinus cavities are included %
			if NSIN > 0 & (j == Nsin(1) | j == Nsin(2))
				if j == Nsin(1)
					k = 1;
				end
				if j == Nsin(2)
					k = 2;
				end
				tmp = Rsin(k) + w*Lsin(k)*i;
				tmp = tmp + 1/(0.0 + w*Csin(k)*i);
				a11 = 1.0 + 0.0*i;
				a12 = 0.0 + 0.0*i;
				a21 = -1/tmp;
				a22 = a11;
				An = a11*Anold + a12*Cnold;
				Bn = a11*Bnold + a12*Dnold;
				Cn = a21*Anold + a22*Cnold;
				Dn = a21*Bnold + a22*Dnold;
				Anold = An;
				Bnold = Bn;
				Cnold = Cn;
				Dnold = Dn;
			end
		end
		Znr = (0.0 + w*RNr*LNr*i)/(RNr + w*LNr*i);
		Hnum = Bn - Dn*Znr;
		Hdem = Cn*Znr - An;
		Zn = Hnum / Hdem;

% from glottis to velopharyngeal port %
		Avpold = 1.0 + 0.0*i;
		Bvpold = 0.0 + 0.0*i;
		Cvpold = 0.0 + 0.0*i;
		Dvpold = 1.0 + 0.0*i;
		for j=1:NTS,
			z = Rs(j) + w*Ls(j)*i;
			wa = Rw(j) + w*Lw(j)*i;
			Zwall = Gp(j) + w*Cp(j)*i;
			y = Zwall + 1/wa;
			Z0 = (z/y)^0.5;
			r = (z*y)^0.5;
			Za = Z0*tanh(leng(j)/2.0*r);
			Zb = Z0/sinh(leng(j)*r);
			a11 = 1.0 + Za/Zb;
			a12 = -(2.0*Za + Za*(Za/Zb));
			a21 = -1/Zb;
			a22 = a11;
			Avp = a11*Avpold + a12*Cvpold;
			Bvp = a11*Bvpold + a12*Dvpold;
			Cvp = a21*Avpold + a22*Cvpold;
			Dvp = a21*Bvpold + a22*Dvpold;
			Avpold = Avp;
			Bvpold = Bvp;
			Cvpold = Cvp;
			Dvpold = Dvp;
		end
		Zno = (Zn*Zo)/(Zn+Zo);
		Hnum = Bvp - Dvp*Zno;
		Hdem = Cvp*Zno - Avp;
		Zp = Hnum / Hdem;

% combine the previous three parts and subglottal system if included %
		a11 = 1.0 + 0.0*i;
		a12 = 0.0 + 0.0*i;
		a21 = -1/Zn;
		a22 = a11;
		Ai = a11*Avp + a12*Cvp;
		Bi = a11*Bvp + a12*Dvp;
		Ci = a21*Avp + a22*Cvp;
		Di = a21*Bvp + a22*Dvp;
		Aiold = Ai;
		Biold = Bi;
		Ciold = Ci;
		Diold = Di;
		Ai = Ao*Aiold + Bo*Ciold;
		Bi = Ao*Biold + Bo*Diold;
		Ci = Co*Aiold + Do*Ciold;
		Di = Co*Biold + Do*Diold;
		Hnum = (Bi*Ci) - (Ai*Di);
		Hdem = (Ci*Zor) - Ai;
		if glotoc == 1
			Hv = Hnum/Hdem;
		else
			Zsubf(l) = Zsubf(l)/(Zp + Zsubf(l));
			Hv = Zsubf(l)*(Hnum/Hdem);
		end

		a11 = 1.0 + 0.0*i;
		a12 = 0.0 + 0.0*i;
		a21 = -1/Zo;
		a22 = a11;
		Ai = a11*Avp + a12*Cvp;
		Bi = a11*Bvp + a12*Dvp;
		Ci = a21*Avp + a22*Cvp;
		Di = a21*Bvp + a22*Dvp;
		Aiold = Ai;
		Biold = Bi;
		Ciold = Ci;
		Diold = Di;
		Ai = An*Aiold + Bn*Ciold;
		Bi = An*Biold + Bn*Diold;
		Ci = Cn*Aiold + Dn*Ciold;
		Di = Cn*Biold + Dn*Diold;
		Hnum = (Bi*Ci) - (Ai*Di);
		Hdem = (Ci*Znr) - Ai;
		if glotoc == 1
			Hn = Hnum / Hdem;
		else
			Hn = Zsubf(l)*(Hnum/Hdem);
		end
		TF(l) = Hn + Hv;

		l = l + 1;
		w = w + df;
end