unvoice.m

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

%************************* Turbulence Noise **********************/
if ((exmode(sfno)==1) | (exmode(sfno)==2))
	for j=1:tnnum,
		tnloc = tnloc + j - 1;
		for jj=1:NING,
			tVv(jj*tnnum+j-1)=tQCv(jj*tnnum+j-1)-Ywv(jj)*tQLwv(jj*tnnum+j-1);
		end
		tVv((NING+1)*tnnum+j-1)=tPv((NING+1)*tnnum+j-1)/(smpfrq*Lr)+tVv((NING+1)*tnnum+j-1);
		Fv(2)=bv(1)*tVv(1*tnnum+j-1)+tQLv(2*tnnum+j-1)-bv(2)*tVv(2*tnnum+j-1);
		for jj=2:NTS-1,
			Fv(jj+1) = bv(jj)*tVv(jj*tnnum+j-1) + tQLv((jj+1)*tnnum+j-1) - bv(jj+1)*tVv((jj+1)*tnnum+j-1);
		end
		tFnc = tQLnts(j) + bv(NTS)*tVv(NTS*tnnum+j-1);
		Fv(NTS+1)=tQLv((NTS+1)*tnnum+j-1)-bv(NTS+1)*tVv((NTS+1)*tnnum+j-1);
		for jj=NTS+2:NING,
			Fv(jj)=tQLv(jj*tnnum+j-1)+bv(jj-1)*tVv((jj-1)*tnnum+j-1)-bv(jj)*tVv(jj*tnnum+j-1);
		end
		Fv(NING+1)=bv(NING)*tVv(NING*tnnum+j-1)+bv(NING+1)*tVv((NING+1)*tnnum+j-1)+tQLning(j);
		Fv(tnloc-1) = bv(tnloc-2)*tVv((tnloc-2)*tnnum+j-1)-bv(tnloc-1)*tVv((tnloc-1)*tnnum+j-1)+tQLv((tnloc-1)*tnnum+j-1)-bv(tnloc-1)*Uvvn(j);
		Fv(tnloc) = bv(tnloc-1)*tVv((tnloc-1)*tnnum+j-1)-bv(tnloc)*tVv(tnloc*tnnum+j-1)+tQLvn(j)+tQLv(tnloc*tnnum+j-1)+(av(tnloc-1)+bv(tnloc-1))*Uvvn(j);
		for jj=1:NTsno,
			tVn(jj*tnnum+j-1)=tQCn(jj*tnnum+j-1)-Ywn(jj)*tQLwn(jj*tnnum+j-1);
		end
		tVn((NTsno+1)*tnnum+j-1) = tPn((NTsno+1)*tnnum+j-1)/(smpfrq*LNr)+tVn((NTsno+1)*tnnum+j-1);
		Fn(1) = tQLn(1*tnnum+j-1) - bn(1)*tVn(1*tnnum+j-1);
		for jj=2:NTsno,
			Fn(jj) = tQLn(jj*tnnum+j-1)+bn(jj-1)*tVn((jj-1)*tnnum+j-1)-bn(jj)*tVn(jj*tnnum+j-1);
		end
		Fn(NTsno+1)=bn(NTsno)*tVn(NTsno*tnnum+j-1)+bn(NTsno+1)*tVn((NTsno+1)*tnnum+j-1)+tQLntsno(j);

	  	if (NSIN == 1) 
			Fsin(1) = tQLsin(1*tnnum+j-1)+tQLns(1*tnnum+j-1)+bn(Nsin(1)+1)*tVn((Nsin(1)+1)*tnnum+j-1)-tVCsin(1*tnnum+j-1);
			Fn(Nsin(1)+2) = tQLn((Nsin(1)+2)*tnnum+j-1)-tQLsin(1*tnnum+j-1)+tVCsin(1*tnnum+j-1)-bn(Nsin(1)+2)*tVn((Nsin(1)+2)*tnnum+j-1);
		end
		if (NSIN == 2) 
			Fsin(1) = tQLsin(1*tnnum+j-1)+tQLns(1*tnnum+j-1)+bn(Nsin(1)+1)*tVn((Nsin(1)+1)*tnnum+j-1)-tVCsin(1*tnnum+j-1);
			Fn(Nsin(1)+2) = tQLn((Nsin(1)+2)*tnnum+j-1)-tQLsin(1*tnnum+j-1)+tVCsin(1*tnnum+j-1)-bn(Nsin(1)+2)*tVn((Nsin(1)+2)*tnnum+j-1);
			Fsin(2) = tQLsin(2*tnnum+j-1)+tQLns(2*tnnum+j-1)+bn(Nsin(2)+1)*tVn((Nsin(2)+1)*tnnum+j-1)-tVCsin(2*tnnum+j-1);
			Fn(Nsin(2)+2)=tQLn((Nsin(2)+2)*tnnum+j-1)-tQLsin(2*tnnum+j-1)+tVCsin(2*tnnum+j-1)-bn(Nsin(2)+2)*tVn((Nsin(2)+2)*tnnum+j-1);
		end

% SOLVE by Gauss elimination */
% ********* Nasal Tract *********/
		if (NSIN==0)
			for jj=1:NTsno+1,
				HNs(jj) = HN(jj);
				bns(jj) = bn(jj);
				Fns(jj) = Fn(jj);
			end
		elseif (NSIN==1)
			for jj=1:Nsin(1)+1,
				HNs(jj) = HN(jj);
				bns(jj) = bn(jj);
				Fns(jj) = Fn(jj);
			end
			HNs(Nsin(1)+2) = Hsf(1);
			bns(Nsin(1)+2) = bsin(1);
			Fns(Nsin(1)+2) = Fsin(1);
			for jj=Nsin(1)+3:NTsno+NSIN+1,
				if (jj > Nsin(1)+3)
					HNs(jj) = HN(jj-1);
				else
					HNs(jj) = Hsb(1);
				end
				bns(jj) = bn(jj-1);
				Fns(jj) = Fn(jj-1);
			end
		elseif (NSIN==2)
			for jj=1:Nsin(1)+1,
				HNs(jj) = HN(jj);
				bns(jj) = bn(jj);
				Fns(jj) = Fn(jj);
			end
			HNs(Nsin(1)+2) = Hsf(1);
			bns(Nsin(1)+2) = bsin(1);
			Fns(Nsin(1)+2) = Fsin(1);
			for jj=Nsin(1)+3:Nsin(2)+2,
				if (jj > Nsin(1)+3)
					HNs(jj) = HN(jj-1);
				else
					HNs(jj) = Hsb(1);
				end
				bns(jj) = bn(jj-1);
				Fns(jj) = Fn(jj-1);
			end
			HNs(Nsin(2)+3) = Hsf(2);
			bns(Nsin(2)+3) = bsin(2);
			Fns(Nsin(2)+3) = Fsin(2);
			for jj=Nsin(2)+4:NTsno+NSIN+1,
				if (jj > Nsin(2)+4)
					HNs(jj) = HN(jj-2);
				else
					HNs(jj) = Hsb(2);
				end
				HNs(jj) = HN(jj-2);
				bns(jj) = bn(jj-2);
				Fns(jj) = Fn(jj-2);
			end
		end
		QQn(NTsno+NSIN+1) = HNs(NTsno+NSIN+1);
		RRn(NTsno+NSIN+1) = Fns(NTsno+NSIN+1);
		for jj=NTsno+NSIN:-1:1,
			QQn(jj) = HNs(jj) - bns(jj)*bns(jj)/QQn(jj+1);
			RRn(jj) = Fns(jj) + bns(jj)*RRn(jj+1)/QQn(jj+1);
		end
%******** Oral Tract *********/
		QQv(NING+1) = HV(NING+1);
		RRv(NING+1) = Fv(NING+1);
		for jj=NING:-1:NTS+1,
			QQv(jj) = HV(jj) - bv(jj)*bv(jj)/QQv(jj+1);
			RRv(jj) = Fv(jj) + bv(jj)*RRv(jj+1)/QQv(jj+1);
		end
%******** Pharyngeal Tract *********/
		QQv(2) = HV(2);
		RRv(2) = Fv(2);
		for jj=3:NTS,
			QQv(jj) = HV(jj) - bv(jj-1)*bv(jj-1)/QQv(jj-1);
			RRv(jj) = Fv(jj) + bv(jj-1)*RRv(jj-1)/QQv(jj-1);
		end
		QNC = tHnc - bv(NTS)*bv(NTS)/QQv(NTS);
		RNC = tFnc + bv(NTS)*RRv(NTS)/QQv(NTS);
%******** Compute Unc and Pnc *************/
		num1 = RNC/QNC - RRv(NTS+1)/QQv(NTS+1) - RRn(1)/QQn(1);
		dem1 = 1.0/QNC + 1.0/QQv(NTS+1) + 1.0/QQn(1);
		Pnc = num1 / dem1;
		Unc = (RNC - Pnc) / QNC;
%******** Backward substitution procedure *********/
% Nasal tract */
		Uns(1) = (RRn(1) + Pnc) / QQn(1);
		for jj=2:NTsno+NSIN+1,
			Uns(jj) = (RRn(jj) + bns(jj-1)*Uns(jj-1)) / QQn(jj);
		end
		if (NSIN==0)
			for jj=1:NTsno+1,
				tUn(jj*tnnum+j-1) = Uns(jj);
			end
		elseif (NSIN==1)
			for jj=1:Nsin(1)+1,
				tUn(jj*tnnum+j-1) = Uns(jj);
			end
			tUsin(1*tnnum+j-1) = Uns(Nsin(1)+2);
			for jj=Nsin(1)+2:NTsno+1,
				tUn(jj*tnnum+j-1) = Uns(jj+1);
			end
		elseif (NSIN==2)
			for jj=1:Nsin(1)+1,
				tUn(jj*tnnum+j-1) = Uns(jj);
			end
			tUsin(1*tnnum+j-1) = Uns(Nsin(1)+2);
			for jj=Nsin(1)+2:Nsin(2)+2,
				tUn(jj*tnnum+j-1) = Uns(jj+1);
			end
			tUsin(2*tnnum+j-1) = Uns(Nsin(2)+3);
			for jj=Nsin(2)+3:NTsno+1,
				tUn(jj*tnnum+j-1) = Uns(jj+2);
			end
		end

% Oral tract */
		tUv((NTS+1)*tnnum+j-1) = (RRv(NTS+1) + Pnc) / QQv(NTS+1);
		for jj=NTS+2:NING+1,
			tUv(jj*tnnum+j-1)=(RRv(jj)+bv(jj-1)*tUv((jj-1)*tnnum+j-1))/QQv(jj);
		end
		tUvp(j) = tUv(tnloc*tnnum+j-1) - Uvvn(j);

% Pharyngeal tract */
		tUv(NTS*tnnum+j-1) = (RRv(NTS) + bv(NTS)*Unc) / QQv(NTS);
		for jj=NTS-1:-1:2,
			tUv(jj*tnnum+j-1) = (RRv(jj) + bv(jj)*tUv((jj+1)*tnnum+j-1))/QQv(jj);
		end
		tPv(1*tnnum+j-1) = bv(1)*( - tUv(2*tnnum+j-1) + tVv(1*tnnum+j-1));
		for jj=2:NTS-1,
			tPv(jj*tnnum+j-1)=bv(jj)*(tUv(jj*tnnum+j-1)-tUv((jj+1)*tnnum+j-1) + tVv(jj*tnnum+j-1));
		end
		tPv(NTS*tnnum+j-1)=bv(NTS)*(tUv(NTS*tnnum+j-1)-Unc+tVv(NTS*tnnum+j-1));
		for jj=NTS+1:NING,
			tPv(jj*tnnum+j-1)=bv(jj)*(tUv(jj*tnnum+j-1)-tUv((jj+1)*tnnum+j-1)+tVv(jj*tnnum+j-1));
		end
		tPv((tnloc-1)*tnnum+j-1)=bv(tnloc-1)*(tUv((tnloc-1)*tnnum+j-1)-tUvp(j) + tVv((tnloc-1)*tnnum+j-1));
		tPv((NING+1)*tnnum+j-1) = bv(NING+1)*(tUv((NING+1)*tnnum+j-1)- tVv((NING+1)*tnnum+j-1));

		for jj=1:NTsno,
			tPn(jj*tnnum+j-1)=bn(jj)*(tUn(jj*tnnum+j-1)-tUn((jj+1)*tnnum+j-1)+ tVn(jj*tnnum+j-1));
		end
		tPn((NTsno+1)*tnnum+j-1) = bn(NTsno+1)*(tUn((NTsno+1)*tnnum+j-1)-tVn((NTsno+1)*tnnum+j-1));
		
		if (NSIN == 1)
			tPn((Nsin(1)+1)*tnnum+j-1)=bn(Nsin(1)+1)*(tUn((Nsin(1)+1)*tnnum+j-1)-tUsin(1*tnnum+j-1) + tVn((Nsin(1)+1)*tnnum+j-1));
			tPsin(1*tnnum+j-1)=bsin(1)*(tUsin(1*tnnum+j-1)-tUn((Nsin(1)+2)*tnnum+j-1))-tQLsin(1*tnnum+j-1)+tVCsin(1*tnnum+j-1);
		end
		if (NSIN == 2) 
			tPn((Nsin(1)+1)*tnnum+j-1) = bn(Nsin(1)+1)*(tUn((Nsin(1)+1)*tnnum+j-1)-tUsin(1*tnnum+j-1)+tVn((Nsin(1)+1)*tnnum+j-1));
			tPsin(1*tnnum+j-1)=bsin(1)*(tUsin(1*tnnum+j-1)-tUn((Nsin(1)+2)*tnnum+j-1))-tQLsin(1*tnnum+j-1)+tVCsin(1*tnnum+j-1);
			tPn((Nsin(2)+1)*tnnum+j-1)=bn(Nsin(2)+1)*(tUn((Nsin(2)+1)*tnnum+j-1)-tUsin(2*tnnum+j-1)+tVn((Nsin(2)+1)*tnnum+j-1));
			tPsin(2*tnnum+j-1)=bsin(2)*(tUsin(2*tnnum+j-1)-tUn((Nsin(2)+2)*tnnum+j-1))-tQLsin(2*tnnum+j-1)+tVCsin(2*tnnum+j-1);
		end


% refresh QLv, QCv, QLwv, u3v, QLn, QCn, QLwn, u3n,..., etc. */

		for jj=2:NTS,
			tQLv(jj*tnnum+j-1)=4.0*smpfrq*(Ls(jj-1)/2.0*leng(jj-1)+Ls(jj)/2.0*leng(jj))*tUv(jj*tnnum+j-1)-tQLv(jj*tnnum+j-1);
		end
		tQLnts(j)=4.0*smpfrq*Ls(NTS)/2.0*leng(NTS)*Unc-tQLnts(j);
		tQLv((NTS+1)*tnnum+j-1)=4.0*smpfrq*Ls(NTS+1)/2.0*leng(NTS+1)*tUv((NTS+1)*tnnum+j-1) - tQLv((NTS+1)*tnnum+j-1);
		for jj=NTS+2:NING,
			if (jj==tnloc)
				tQLv(tnloc*tnnum+j-1) = 4.0*smpfrq*Ls(tnloc)/2.0*leng(tnloc)*tUv(tnloc*tnnum+j-1)-tQLv(tnloc*tnnum+j-1);
			else
				tQLv(jj*tnnum+j-1) = 4.0*smpfrq*(Ls(jj-1)/2.0*leng(jj-1) +Ls(jj)/2.0*leng(jj))*tUv(jj*tnnum+j-1)-tQLv(jj*tnnum+j-1);
			end
		end

		tQLvn(j) = 4.0*smpfrq*Ls(tnloc-1)/2.0*leng(tnloc-1)*tUvp(j) - tQLvn(j);
		tQLning(j) = 4.0*smpfrq*Ls(NING)/2.0*leng(NING)*tUv((NING+1)*tnnum+j-1) - tQLning(j);

		for jj=1:NING,
			tQCv(jj*tnnum+j-1)=4.0*smpfrq*Cp(jj)*leng(jj)*tPv(jj*tnnum+j-1)-tQCv(jj*tnnum+j-1);
			u3v(jj) = Ywv(jj)*(tPv(jj*tnnum+j-1) + tQLwv(jj*tnnum+j-1));
			tQLwv(jj*tnnum+j-1)=4.0*smpfrq*Lw(jj)*u3v(jj)-tQLwv(jj*tnnum+j-1);
		end

		tQLn(1*tnnum+j-1) = 4.0*smpfrq*LNs(1)/2.0*NTslen*tUn(1*tnnum+j-1)-tQLn(1*tnnum+j-1);
		for jj=2:NTsno,
			tQLn(jj*tnnum+j-1) = 4.0*smpfrq*(LNs(jj-1)/2.0+LNs(jj)/2.0)*NTslen*tUn(jj*tnnum+j-1) - tQLn(jj*tnnum+j-1);
		end
		tQLntsno(j) = 4.0*smpfrq*LNs(NTsno-1)/2.0*NTslen*tUn((NTsno+1)*tnnum+j-1)- tQLntsno(j);

		for jj=1:NTsno,
			tQCn(jj*tnnum+j-1) = 4.0*smpfrq*CNp(jj)*NTslen*tPn(jj*tnnum+j-1) - tQCn(jj*tnnum+j-1);
			u3n(jj) = Ywn(jj)*(tPn(jj*tnnum+j-1) + tQLwn(jj*tnnum+j-1));
			tQLwn(jj*tnnum+j-1) = 4.0*smpfrq*LNw(jj)*u3n(jj) - tQLwn(jj*tnnum+j-1);
		end
		if (NSIN == 1) 
			u3sin(1)=(tPsin(1*tnnum+j-1)+tQLsin(1*tnnum+j-1)-tVCsin(1*tnnum+j-1))/bsin(1);
			tQLsin(1*tnnum+j-1) = 4.0*smpfrq*Lsin(1)*u3sin(1)-tQLsin(1*tnnum+j-1);
			tVCsin(1*tnnum+j-1) = u3sin(1)/(smpfrq*Csin(1))+ tVCsin(1*tnnum+j-1);
			tQLns(1*tnnum+j-1) = 4.0*smpfrq*LNs(Nsin(1))/2.0*NTslen*tUsin(1*tnnum+j-1)-tQLns(1*tnnum+j-1);
			tQLn((Nsin(1)+2)*tnnum+j-1) = 4.0*smpfrq*LNs(Nsin(1)+1)/2.0*NTslen*tUn((Nsin(1)+2)*tnnum+j-1)- tQLn((Nsin(1)+2)*tnnum+j-1);
		end
		if (NSIN == 2) 
			u3sin(1) = (tPsin(1*tnnum+j-1)+tQLsin(1*tnnum+j-1)-tVCsin(1*tnnum+j-1))/bsin(1);
			tQLsin(1*tnnum+j-1) = 4.0*smpfrq*Lsin(1)*u3sin(1)-tQLsin(1*tnnum+j-1);
			tVCsin(1*tnnum+j-1) = u3sin(1)/(smpfrq*Csin(1))+ tVCsin(1*tnnum+j-1);
			tQLns(1*tnnum+j-1) = 4.0*smpfrq*LNs(Nsin(1))/2.0*NTslen*tUsin(1*tnnum+j-1)-tQLns(1*tnnum+j-1);
			tQLn((Nsin(1)+2)*tnnum+j-1) = 4.0*smpfrq*LNs(Nsin(1)+1)/2.0*NTslen*tUn((Nsin(1)+2)*tnnum+j-1)- tQLn((Nsin(1)+2)*tnnum+j-1);
			u3sin(2) = (tPsin(2*tnnum+j-1)+tQLsin(2*tnnum+j-1)-tVCsin(2*tnnum+j-1))/bsin(2);
			tQLsin(2*tnnum+j-1) = 4.0*smpfrq*Lsin(2)*u3sin(2)-tQLsin(2*tnnum+j-1);
			tVCsin(2*tnnum+j-1) = u3sin(2)/(smpfrq*Csin(2))+ tVCsin(2*tnnum+j-1);
			tQLns(2*tnnum+j-1) = 4.0*smpfrq*LNs(Nsin(2))/2.0*NTslen*tUsin(2*tnnum+j-1)-tQLns(2*tnnum+j-1);
			tQLn((Nsin(2)+2)*tnnum+j-1) = 4.0*smpfrq*LNs(Nsin(2)+1)/2.0*NTslen*tUn((Nsin(2)+2)*tnnum+j-1) - tQLn((Nsin(2)+2)*tnnum+j-1);
		end

	end
end