⭐ 欢迎来到虫虫下载站! | 📦 资源下载 📁 资源专辑 ℹ️ 关于我们
⭐ 虫虫下载站

📄 syncaltf.m

📁 这是一个用于语音信号处理的工具箱
💻 M
📖 第 1 页 / 共 2 页
字号:
% syncaltf.m

clear i;
	glotoc = 1;
	if exist('EX_loc')
		exloc=EX_loc(vfno);
	else
		exloc = 0;
	end
	if (jsaspsg(sfno)==4 | jsaspsg(sfno)==5 | jsaspsg(sfno)==6 | jsaspsg(sfno)==7)
		subglimp;
		glotoc=0;
	end

	NING = 60;
	ISEC = 60/NING;
	NTS = round(29/ISEC);
	Nfft = 64;
	
	NTlen = 11;
	NTsno = 11;
	NTslen = NTlen/NTsno;
	
	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;
	NTArea(1) = nt1atemp + 0.0001;
	del_n = (NTArea(4) - NTArea(1))/3.0;
	NTArea(2) = NTArea(1) + del_n;
	NTArea(3) = NTArea(2) + del_n;

%	NTArea(2) = 2.0;
%	NTArea(3) = 3.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,
 
%   Perimeter = 2.0 * sqrt(pi * A) <----- circle
%            = 4.0 * sqrt(pi * A) <----- ellipse

		peri = 4.0 * sqrt(pi*areatemp(j));

		Rs(j) = peri * rdem/((areatemp(j))^2);
		Ls(j) = rho / areatemp(j);
		Cp(j) = areatemp(j) / rc2;
		Gp(j) = peri * Gpc;

		wall = sqrt(pi)*areatemp(j)*lengtemp(j);        % Maeda */
%               wall = peri*sl->dl[i];  Flanagan */
		Rw(j) = viscous / wall;
		Lw(j) = mass / wall;
%               Cw[i] = wall / compilance;      */
	end
	Rr = (128.0*rho*c)/(9.0*pi*pi*areatemp(NING));
	Lr = 1.5 * ( 8.0*rho / (3.0*pi*sqrt(pi*areatemp(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,
 
%   Perimeter = 2.0 * sqrt(pi * A) <----- circle
%            = 4.0 * sqrt(pi * A) <----- ellipse

		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 */
%               wall = peri*NTslen; Flanagan */
		RNw(j) = viscous / wall;
		LNw(j) = mass / wall;
%               CNw(j) = wall / compilance;     */
	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.

% function TFcaseA(tfra, glotoc, exloc, Nfft)

if exloc == 0

% When the excitation located at the glottis.
%	TFcaseA(tfra, glotoc, exloc, Nfft);	

	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
		tmp = Hn + Hv;
		TF(l) = 20.0*log10(abs(tmp));
%		TF(l) = Hn + Hv;

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

	elseif exloc > 0 & exloc < NTS

% When the excitation located at the pharyngeal tract
%	TFcaseB(tfra, glotoc, exloc, Nfft);

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

	while l <= Nfft,

% from glottis to excitation location. Also considered subglottal system %
		Avpold = 1 + 0*i;
		Bvpold = 0 + 0*i;
		Cvpold = 0 + 0*i;
		Dvpold = 1 + 0*i;
		for j=1:exloc,
			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
		if glotoc == 1
			Zsubf(l) = Bvp/Dvp;
		else
			Hnum = Avp*Zsubf(l) + Bvp;
			Hdem = Cvp*Zsubf(l) + Dvp;
			Zsubf(l) = Hnum/Hdem;
		end

% from excitation location 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=exloc+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

% 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 + 0.0*i;
		Bnold = 0.0 + 0.0*i;
		Cnold = 0.0 + 0.0*i;
		Dnold = 1.0 + 0.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 for sinus cavities coupling %
			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;

⌨️ 快捷键说明

复制代码 Ctrl + C
搜索代码 Ctrl + F
全屏模式 F11
切换主题 Ctrl + Shift + D
显示快捷键 ?
增大字号 Ctrl + =
减小字号 Ctrl + -