syntest1.m
来自「这是一个用于语音信号处理的工具箱」· M 代码 · 共 526 行
M
526 行
% syntest1.m
% modified by D. G. Childers 8/18/98
ient = 1;
% First time called, it needs to initialize force constants (see Maeda, 1982),
% allocate memory, and calculate the impedance related constants for the
% nasal tract from fourth section to nostrils. */
if (ient == 0)
gQLnts = 0.0;
gQLning = 0.0;
gQLntsno = 0.0;
gHnc = 0.0;
gFnc = 0.0;
tHnc = 0.0;
tFnc = 0.0;
for j=1:15,
HNs(j) = 0.0;
bns(j) = 0.0;
Fns(j) = 0.0;
Uns(j) = 0.0;
end
clear QQv;
clear RRv;
clear QQn;
clear RRn;
QQv = zeros(1,NING+2);
RRv = zeros(1,NING+2);
QQn = zeros(1,NTsno+NSIN+2);
RRn = zeros(1,NTsno+NSIN+2);
clear av;
clear bv;
clear Ywv;
clear Fv;
clear HV;
clear u3v;
clear gUv;
clear gPv;
clear gQLv;
clear gVv;
clear gQCv;
clear gQLwv;
gQCv=zeros(1,60);
gQLwv=zeros(1,60);
gPv=zeros(1,61);
if (exist('tUv') & tnnum >= 1)
clear tUv;
clear tPv;
clear tQLv;
clear tVv;
clear tQCv;
clear tQLwv;
end
av = zeros(1,NING+1);
bv = zeros(1,NING+2);
Ywv = zeros(1,NING+1);
Fv = zeros(1,NING+2);
HV = zeros(1,NING+2);
u3v = zeros(1,NING+1);
gPv = zeros(1,NING+2);
gUv = zeros(1,NING+2);
gQLv = zeros(1,NING+1);
gVv = zeros(1,NING+2);
gQCv = zeros(1,NING+1);
gQLwv= zeros(1,NING+1);
tPv = zeros(1,NING+2);
tUv = zeros(1,NING+2);
tQLv = zeros(1,NING+1);
tVv = zeros(1,NING+2);
tQCv = zeros(1,NING+1);
tQLwv= zeros(1,NING+1);
clear an;
clear bn;
clear Ywn;
clear Fn;
clear HN;
clear u3n;
clear gUn;
clear gPn;
clear gQLn;
clear gVn;
clear gQCn;
clear gQLwn;
if (exist('tUn') & tnnum >= 1)
clear tUn;
clear tPn;
clear tQLn;
clear tVn;
clear tQCn;
clear tQLwn;
end
an = zeros(1,NTsno+1);
bn = zeros(1,NTsno+2);
Ywn = zeros(1,NTsno+1);
Fn = zeros(1,NTsno+2);
HN = zeros(1,NTsno+2);
u3n = zeros(1,NTsno+1);
gUn = zeros(1,NTsno+2);
gPn = zeros(1,NTsno+2);
gQLn = zeros(1,NTsno+1);
gVn = zeros(1,NTsno+2);
gQCn = zeros(1,NTsno+1);
gQLwn= zeros(1,NTsno+1);
tUn = zeros(1,(NTsno+2)*tnnum);
tPn = zeros(1,(NTsno+2)*tnnum);
tQLn = zeros(1,(NTsno+1)*tnnum);
tVn = zeros(1,(NTsno+2)*tnnum);
tQCn = zeros(1,(NTsno+1)*tnnum);
tQLwn= zeros(1,(NTsno+1)*tnnum);
clear tQLvn;
clear tUvp;
clear tQLnts;
clear tQLning;
clear tQLntsno;
tQLvn = zeros(1,tnnum+1);
tUvp = zeros(1,tnnum+1);
tQLnts = zeros(1,tnnum+1);
tQLning = zeros(1,tnnum+1);
tQLntsno = zeros(1,tnnum+1);
if (NSIN > 0)
clear bsin;
clear Hsb;
clear Hsf;
clear Fsin;
clear u3sin;
clear gQLsin;
clear gVCsin;
clear gUsin;
clear gPsin;
clear gQLns;
if (exist('tQLsin') & tnnum >= 1)
clear tQLsin;
clear tVCsin;
clear tUsin;
clear tPsin;
clear tQLns;
end
bsin = zeros(1,NSIN+1);
Hsb = zeros(1,NSIN+1);
Hsf = zeros(1,NSIN+1);
Fsin = zeros(1,NSIN+1);
u3sin = zeros(1,NSIN+1);
gQLsin = zeros(1,NSIN+1);
gVCsin = zeros(1,NSIN+1);
gUsin = zeros(1,NSIN+1);
gPsin = zeros(1,NSIN+1);
gQLns = zeros(1,NSIN+1);
tQLsin=zeros(1,(NSIN+1)*tnnum);
tVCsin=zeros(1,(NSIN+1)*tnnum);
tUsin=zeros(1,(NSIN+1)*tnnum);
tPsin=zeros(1,(NSIN+1)*tnnum);
tQLns=zeros(1,(NSIN+1)*tnnum);
end
for j=4:NTsno,
an(j) = 0.5*(2.0*smpfrq*LNs(j) + RNs(j))*NTslen;
Ywn(j) = 1.0/((2.0*smpfrq*LNw(j) + RNw(j)));
bn(j) = 1.0 / ((2.0*smpfrq*CNp(j)+GNp(j))*NTslen + Ywn(j));
if (j > 4)
HN(j) = an(j-1) + bn(j-1) + an(j) + bn(j);
end
end
bn(NTsno+1) = 1.0 /(1.0 / RNr + 1.0 / (2.0*smpfrq*LNr));
HN(NTsno+1) = an(NTsno) + bn(NTsno) + bn(NTsno+1);
end %% end of ient = 0 */
%% Calculate the elements of acoustic matrix equations */
if (ient < 2)
for j=1:NING,
av(j) = 0.5*(2.0*smpfrq*Ls(j) + Rs(j))*lengtemp(j);
Ywv(j) = 1.0 / ((2.0*smpfrq*Lw(j) + Rw(j)));
bv(j) = 1.0 / ((2.0*smpfrq*Cp(j) + Gp(j))*lengtemp(j) + Ywv(j));
if (j==1 | j==NTS+1)
HV(j) = av(j) + bv(j);
else
HV(j) = av(j-1) + bv(j-1) + av(j) + bv(j);
end
end
gHnc = av(NTS) + bv(NTS);
tHnc = av(NTS) + bv(NTS);
bv(NING+1) = 1.0 /(1.0 / Rr + 1.0 / (2.0*smpfrq*Lr));
HV(NING+1) = av(NING) + bv(NING) + bv(NING+1);
for j=1:3,
an(j) = 0.5*(2.0*smpfrq*LNs(j) + RNs(j))*NTslen;
Ywn(j) = 1.0/((2.0*smpfrq*LNw(j) + RNw(j)));
bn(j) = 1.0 / ((2.0*smpfrq*CNp(j)+GNp(j))*NTslen + Ywn(j));
if (j > 1)
HN(j) = an(j-1) + bn(j-1) + an(j) + bn(j);
else
HN(j) = an(j) + bn(j);
end
end
HN(4) = an(3) + bn(3) + an(4) + bn(4);
if (NSIN == 1)
bsin(1)=2.0*smpfrq*Lsin(1)+Rsin(1)+0.5/(smpfrq*Csin(1));
Hsf(1) = an(Nsin(1)+1) + bn(Nsin(1)+1) + bsin(1);
Hsb(1) = an(Nsin(1)+2) + bn(Nsin(1)+2) + bsin(1);
end
if (NSIN == 2)
bsin(1)=2.0*smpfrq*Lsin(1)+Rsin(1)+0.5/(smpfrq*Csin(1));
Hsf(1) = an(Nsin(1)+1) + bn(Nsin(1)+1) + bsin(1);
Hsb(1) = an(Nsin(1)+2) + bn(Nsin(1)+2) + bsin(1);
bsin(2)=2.0*smpfrq*Lsin(2)+Rsin(2)+0.5/(smpfrq*Csin(2));
Hsf(2) = an(Nsin(2)+1) + bn(Nsin(2)+1) + bsin(2);
Hsb(2) = an(Nsin(2)+2) + bn(Nsin(2)+2) + bsin(2);
end
end
pulsesamp_len=length(pulsesamp);
if (exmode(sfno)==1) % turbulence source
gUv(1) = TN_vvdc(sfno) + 200.0*(rand(1)-0.5);
elseif (exmode(sfno)==2 & srcidx<=pulsesamp_len) % mixed
gUv(1) = pulsesamp(srcidx) + TN_vvdc(sfno) + 200.0*(rand(1)-0.5);
srcidx = srcidx + 1;
elseif srcidx<=pulsesamp_len % glottal source
gUv(1) = pulsesamp(srcidx);
srcidx = srcidx + 1;
end
% if (exmode(sfno)==1) % turbulence source
% gUv(1) = TN_vvdc(sfno) + 200.0*(rand(1)-0.5);
% elseif (exmode(sfno)==2) % mixed
% gUv(1) = pulsesamp(srcidx) + TN_vvdc(sfno) + 200.0*(rand(1)-0.5);
% srcidx = srcidx + 1;
% else % glottal source
% gUv(1) = pulsesamp(srcidx);
% srcidx = srcidx + 1;
% end
for j=1:NING,
gVv(j) = gQCv(j) - Ywv(j)*gQLwv(j);
end
gVv(NING+1) = gPv(NING+1) / (smpfrq*Lr) + gVv(NING+1);
%% Calculate the force constants */
Fv(1) = gQLv(1) - bv(1)*gVv(1) - HV(1)*gUv(1);
Fv(2) = bv(1)*gVv(1) + gQLv(2) - bv(2)*gVv(2) + bv(1)*gUv(1);
for j=2:NTS-1,
Fv(j+1) = bv(j)*gVv(j) + gQLv(j+1) - bv(j+1)*gVv(j+1);
end
gFnc = gQLnts + bv(NTS)*gVv(NTS);
Fv(NTS+1) = gQLv(NTS+1) - bv(NTS+1)*gVv(NTS+1);
for j=NTS+2:NING,
Fv(j) = gQLv(j) + bv(j-1)*gVv(j-1) - bv(j)*gVv(j);
end
Fv(NING+1) = bv(NING)*gVv(NING) + bv(NING+1)*gVv(NING+1) + gQLning;
for j=1:NTsno,
gVn(j) = gQCn(j) - Ywn(j)*gQLwn(j);
end
gVn(NTsno+1) = gPn(NTsno+1) / (smpfrq*LNr) + gVn(NTsno+1);
Fn(1) = gQLn(1) - bn(1)*gVn(1);
for j=2:NTsno,
Fn(j) = gQLn(j) + bn(j-1)*gVn(j-1) - bn(j)*gVn(j);
end
Fn(NTsno+1) = bn(NTsno)*gVn(NTsno)+bn(NTsno+1)*gVn(NTsno+1)+gQLntsno;
if (NSIN == 1)
Fsin(1) = gQLsin(1)+gQLns(1)+bn(Nsin(1)+1)*gVn(Nsin(1)+1) - gVCsin(1);
Fn(Nsin(1)+2) = gQLn(Nsin(1)+2)-gQLsin(1)+gVCsin(1) - bn(Nsin(1)+2)*gVn(Nsin(1)+2);
end
if (NSIN == 2)
Fsin(1) = gQLsin(1)+gQLns(1)+bn(Nsin(1)+1)*gVn(Nsin(1)+1) - gVCsin(1);
Fn(Nsin(1)+2) = gQLn(Nsin(1)+2)-gQLsin(1)+gVCsin(1) - bn(Nsin(1)+2)*gVn(Nsin(1)+2);
Fsin(2) = gQLsin(2)+gQLns(2)+bn(Nsin(2)+1)*gVn(Nsin(2)+1) - gVCsin(2);
Fn(Nsin(2)+2) = gQLn(Nsin(2)+2)-gQLsin(2)+gVCsin(2) - bn(Nsin(2)+2)*gVn(Nsin(2)+2);
end
% SOLVE by Gauss elimination*/
%******** Nasal Tract *********/
if (NSIN==0)
% case 0:
for j=1:NTsno+1,
HNs(j) = HN(j);
bns(j) = bn(j);
Fns(j) = Fn(j);
end
% case 1:
elseif (NSIN==1)
for j=1:Nsin(1)+1,
HNs(j) = HN(j);
bns(j) = bn(j);
Fns(j) = Fn(j);
end
HNs(Nsin(1)+2) = Hsf(1);
bns(Nsin(1)+2) = bsin(1);
Fns(Nsin(1)+2) = Fsin(1);
for j=Nsin(1)+3:NTsno+NSIN+1,
if (j > Nsin(1)+3)
HNs(j) = HN(j-1);
else
HNs(j) = Hsb(1);
end
bns(j) = bn(j-1);
Fns(j) = Fn(j-1);
end
% case 2:
elseif (NSIN==2)
for j=1:Nsin(1)+1,
HNs(j) = HN(j);
bns(j) = bn(j);
Fns(j) = Fn(j);
end
HNs(Nsin(1)+2) = Hsf(1);
bns(Nsin(1)+2) = bsin(1);
Fns(Nsin(1)+2) = Fsin(1);
for j=Nsin(1)+3:Nsin(2)+2,
if (j > Nsin(1)+3)
HNs(j) = HN(j-1);
else
HNs(j) = Hsb(1);
end
bns(j) = bn(j-1);
Fns(j) = Fn(j-1);
end
HNs(Nsin(2)+3) = Hsf(2);
bns(Nsin(2)+3) = bsin(2);
Fns(Nsin(2)+3) = Fsin(2);
for j=Nsin(2)+4:NTsno+NSIN+1,
if (j > Nsin(2)+4)
HNs(j) = HN(j-2);
else
HNs(j) = Hsb(2);
end
% HNs(j) = HN(j-2); */
bns(j) = bn(j-2);
Fns(j) = Fn(j-2);
end
end
QQn(NTsno+NSIN+1) = HNs(NTsno+NSIN+1);
RRn(NTsno+NSIN+1) = Fns(NTsno+NSIN+1);
for j=NTsno+NSIN:-1:1,
QQn(j) = HNs(j) - bns(j)*bns(j)/QQn(j+1);
RRn(j) = Fns(j) + bns(j)*RRn(j+1)/QQn(j+1);
end
%******** Oral Tract *********/
QQv(NING+1) = HV(NING+1);
RRv(NING+1) = Fv(NING+1);
for j=NING:-1:NTS+1,
QQv(j) = HV(j) - bv(j)*bv(j)/QQv(j+1);
RRv(j) = Fv(j) + bv(j)*RRv(j+1)/QQv(j+1);
end
%******** Pharyngeal Tract *********/
QQv(1) = -1.0;
RRv(1) = Fv(1);
QQv(2) = HV(2);
RRv(2) = Fv(2);
for j=3:NTS,
QQv(j) = HV(j) - bv(j-1)*bv(j-1)/QQv(j-1);
RRv(j) = Fv(j) + bv(j-1)*RRv(j-1)/QQv(j-1);
end
QNC = gHnc - bv(NTS)*bv(NTS)/QQv(NTS);
RNC = gFnc + 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 j=2:NTsno+NSIN+1,
Uns(j) = (RRn(j) + bns(j-1)*Uns(j-1)) / QQn(j);
end
if (NSIN==0)
% case 0:
for j=1:NTsno+1,
gUn(j) = Uns(j);
end
elseif (NSIN==1)
% case 1:
for j=1:Nsin(1)+1,
gUn(j) = Uns(j);
end
gUsin(1) = Uns(Nsin(1)+2);
for j=Nsin(1)+2:NTsno+1,
gUn(j) = Uns(j+1);
end
% case 2:
elseif (NSIN==2)
for j=1:Nsin(1)+1,
gUn(j) = Uns(j);
end
gUsin(1) = Uns(Nsin(1)+2);
for j=Nsin(1)+2:Nsin(2)+2,
gUn(j) = Uns(j+1);
end
gUsin(2) = Uns(Nsin(2)+3);
for j=Nsin(2)+3:NTsno+1,
gUn(j) = Uns(j+2);
end
end
% Oral tract */
gUv(NTS+1) = (RRv(NTS+1) + Pnc) / QQv(NTS+1);
for j=NTS+2:NING+1,
gUv(j) = (RRv(j) + bv(j-1)*gUv(j-1)) / QQv(j);
end
% Pharyngeal tract */
gUv(NTS) = (RRv(NTS) + bv(NTS)*Unc) / QQv(NTS);
for j=NTS-1:-1:2,
gUv(j) = (RRv(j) + bv(j)*gUv(j+1)) / QQv(j);
end
PE = (RRv(1) + bv(1)*gUv(2)) / QQv(1);
for j=1:NTS-1,
gPv(j) = bv(j)*(gUv(j) - gUv(j+1) + gVv(j));
end
gPv(NTS) = bv(NTS)*(gUv(NTS) - Unc + gVv(NTS));
for j=NTS+1:NING,
gPv(j) = bv(j)*(gUv(j) - gUv(j+1) + gVv(j));
end
gPv(NING+1) = bv(NING+1)*(gUv(NING+1) - gVv(NING+1));
for j=1:NTsno,
gPn(j) = bn(j)*(gUn(j) - gUn(j+1) + gVn(j));
end
gPn(NTsno+1) = bn(NTsno+1)*(gUn(NTsno+1) - gVn(NTsno+1));
if (NSIN == 1)
gPn(Nsin(1)+1) = bn(Nsin(1)+1) * (gUn(Nsin(1)+1)-gUsin(1)+gVn(Nsin(1)+1));
gPsin(1)=bsin(1)*(gUsin(1)-gUn(Nsin(1)+2))-gQLsin(1)+gVCsin(1);
end
if (NSIN == 2)
gPn(Nsin(1)+1) = bn(Nsin(1)+1) * (gUn(Nsin(1)+1)-gUsin(1)+gVn(Nsin(1)+1));
gPsin(1)=bsin(1)*(gUsin(1)-gUn(Nsin(1)+2))-gQLsin(1)+gVCsin(1);
gPn(Nsin(2)+1) = bn(Nsin(2)+1) * (gUn(Nsin(2)+1)-gUsin(2)+gVn(Nsin(2)+1));
gPsin(2)=bsin(2)*(gUsin(2)-gUn(Nsin(2)+2))-gQLsin(2)+gVCsin(2);
end
%% refresh QLv, QCv, QLwv, u3v, QLn, QCn, QLwn, u3n,..., etc. */
gQLv(1) = 4.0*smpfrq*Ls(1)/2.0*lengtemp(1)*gUv(1) - gQLv(1);
for j=2:NTS,
gQLv(j) = 4.0*smpfrq*(Ls(j-1)/2.0*lengtemp(j-1) + Ls(j)/2.0*lengtemp(j))*gUv(j) - gQLv(j);
end
gQLnts = 4.0*smpfrq*Ls(NTS)/2.0*lengtemp(NTS)*Unc - gQLnts;
gQLv(NTS+1) = 4.0*smpfrq*Ls(NTS+1)/2.0*lengtemp(NTS+1)*gUv(NTS+1) - gQLv(NTS+1);
for j=NTS+2:NING,
gQLv(j) = 4.0*smpfrq*(Ls(j-1)/2.0*lengtemp(j-1) + Ls(j)/2.0*lengtemp(j))*gUv(j) - gQLv(j);
end
gQLning = 4.0*smpfrq*Ls(NING)/2.0*lengtemp(NING)*gUv(NING+1) - gQLning;
for j=1:NING,
gQCv(j)=4.0*smpfrq*Cp(j)*lengtemp(j)*gPv(j)-gQCv(j);
u3v(j) = Ywv(j)*(gPv(j) + gQLwv(j));
gQLwv(j) = 4.0*smpfrq*Lw(j)*u3v(j) - gQLwv(j);
end
gQLn(1) = 4.0*smpfrq*LNs(1)/2.0*NTslen*gUn(1) - gQLn(1);
for j=2:NTsno,
gQLn(j) = 4.0*smpfrq*(LNs(j-1)/2.0+LNs(j)/2.0)*NTslen*gUn(j) - gQLn(j);
end
gQLntsno = 4.0*smpfrq*LNs(NTsno)/2.0*NTslen*gUn(NTsno+1) - gQLntsno;
for j=1:NTsno,
gQCn(j) = 4.0*smpfrq*CNp(j)*NTslen*gPn(j) - gQCn(j);
u3n(j) = Ywn(j)*(gPn(j) + gQLwn(j));
gQLwn(j) = 4.0*smpfrq*LNw(j)*u3n(j) - gQLwn(j);
end
if (NSIN == 1)
u3sin(1) = (gPsin(1)+gQLsin(1)-gVCsin(1))/bsin(1);
gQLsin(1) = 4.0*smpfrq*Lsin(1)*u3sin(1) - gQLsin(1);
gVCsin(1) = u3sin(1)/(smpfrq*Csin(1)) + gVCsin(1);
gQLns(1) = 4.0*smpfrq*LNs(Nsin(1))/2.0*NTslen*gUsin(1)-gQLns(1);
gQLn(Nsin(1)+2) = 4.0*smpfrq*LNs(Nsin(1)+1)/2.0*NTslen*gUn(Nsin(1)+2) - gQLn(Nsin(1)+2);
end
if (NSIN == 2)
u3sin(1) = (gPsin(1)+gQLsin(1)-gVCsin(1))/bsin(1);
gQLsin(1) = 4.0*smpfrq*Lsin(1)*u3sin(1) - gQLsin(1);
gVCsin(1) = u3sin(1)/(smpfrq*Csin(1)) + gVCsin(1);
gQLns(1) = 4.0*smpfrq*LNs(Nsin(1))/2.0*NTslen*gUsin(1)-gQLns(1);
gQLn(Nsin(1)+2) = 4.0*smpfrq*LNs(Nsin(1)+1)/2.0*NTslen*gUn(Nsin(1)+2) - gQLn(Nsin(1)+2);
u3sin(2) = (gPsin(2)+gQLsin(2)-gVCsin(2))/bsin(2);
gQLsin(2) = 4.0*smpfrq*Lsin(2)*u3sin(2) - gQLsin(2);
gVCsin(2) = u3sin(2)/(smpfrq*Csin(2)) + gVCsin(2);
gQLns(2) = 4.0*smpfrq*LNs(Nsin(2))/2.0*NTslen*gUsin(2)-gQLns(2);
gQLn(Nsin(2)+2) = 4.0*smpfrq*LNs(Nsin(2)+1)/2.0*NTslen*gUn(Nsin(2)+2) - gQLn(Nsin(2)+2);
end
unvoice;
%% Combine the glottal source with the turbulence source (if exists) */
for j1=1:NING+1,
Uv(j1) = gUv(j1);
Pv(j1) = gPv(j1);
for j=1:tnnum,
Uv(j1) = Uv(j1) + tUv(j1*tnnum+j);
Pv(j1) = Pv(j1) + tPv(j1*tnnum+j);
end
end
for j1=1:NTsno+1,
Un(j1) = gUn(j1);
Pn(j1) = gPn(j1);
for j=1:tnnum,
Un(j1) = Un(j1) + tUn(j1*tnnum+j);
Pn(j1) = Pn(j1) + tPn(j1*tnnum+j);
end
end
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?