synrelo.m
来自「这是一个用于语音信号处理的工具箱」· M 代码 · 共 432 行
M
432 行
% synrelo.m
Unc=0.0; Pnc=0.0; PE=0.0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ient = 0; NSIN = 0;NTsno = 11;srcidx = 1;
% smpfrq = 20000;TNvvdc = 1000;
tnnum = 1;
%exmode = 0;
%tnloc = 53;
%Uvvn(1) = 3;
%Uvvn(2) = 9;
% Excitation at glottis without subglottal system. The glottal source (with
% prefix g) and turbulence source (with prefix t) are solved separately.
%void exgltnos(ient)
gQLnts=0.0; gQLning=0.0; gQLntsno=0.0;
gHnc=0.0; gFnc=0.0;
tHnc=0.0; tFnc=0.0;
% 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;
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)*tnnum);
tUv = zeros(1,(NING+2)*tnnum);
tQLv = zeros(1,(NING+1)*tnnum);
tVv = zeros(1,(NING+2)*tnnum);
tQCv = zeros(1,(NING+1)*tnnum);
tQLwv= zeros(1,(NING+1)*tnnum);
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))*leng(j);
Ywv(j) = 1.0 / ((2.0*smpfrq*Lw(j) + Rw(j)));
bv(j) = 1.0 / ((2.0*smpfrq*Cp(j) + Gp(j))*leng(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(0)+Rsin(0)+0.5/(smpfrq*Csin(0));
Hsf(1) = an(Nsin(0)+1) + bn(Nsin(0)+1) + bsin(1);
Hsb(1) = an(Nsin(0)+2) + bn(Nsin(0)+2) + bsin(1);
end
if (NSIN == 2)
bsin(1)=2.0*smpfrq*Lsin(0)+Rsin(0)+0.5/(smpfrq*Csin(0));
Hsf(1) = an(Nsin(0)+1) + bn(Nsin(0)+1) + bsin(1);
Hsb(1) = an(Nsin(0)+2) + bn(Nsin(0)+2) + bsin(1);
bsin(2)=2.0*smpfrq*Lsin(1)+Rsin(1)+0.5/(smpfrq*Csin(1));
Hsf(2) = an(Nsin(1)+1) + bn(Nsin(1)+1) + bsin(2);
Hsb(2) = an(Nsin(1)+2) + bn(Nsin(1)+2) + bsin(2);
end
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) = 0;
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(0)+1)*gVn(Nsin(0)+1) - gVCsin(1);
Fn(Nsin(0)+2) = gQLn(Nsin(0)+2)-gQLsin(1)+gVCsin(1) - bn(Nsin(0)+2)*gVn(Nsin(0)+2);
end
if (NSIN == 2)
Fsin(1) = gQLsin(1)+gQLns(1)+bn(Nsin(0)+1)*gVn(Nsin(0)+1) - gVCsin(1);
Fn(Nsin(0)+2) = gQLn(Nsin(0)+2)-gQLsin(1)+gVCsin(1) - bn(Nsin(0)+2)*gVn(Nsin(0)+2);
Fsin(2) = gQLsin(2)+gQLns(2)+bn(Nsin(1)+1)*gVn(Nsin(1)+1) - gVCsin(2);
Fn(Nsin(1)+2) = gQLn(Nsin(1)+2)-gQLsin(2)+gVCsin(2) - bn(Nsin(1)+2)*gVn(Nsin(1)+2);
end
F=[];
for j=1:29,
F=[F,Fv(j)];
end
F=[F,gFnc];
for j=30:61,
F=[F,Fv(j)];
end
for j=1:12,
F=[F,Fn(j)];
end
F(75)=0;
F=F';
FFF=zeros(75,75);
FFF(1,1)=-1;
for j=1:29,
FFF(j,j+1)=-bv(j);
end
for j=2:29,
FFF(j,j)=HV(j);
FFF(j+1,j)=-bv(j);
end
FFF(30,30)=gHnc;
FFF(30,31)=1;
FFF(31,31)=-1;
for j=31:61,
FFF(j,j+1)=HV(j-1);
FFF(j,j+2)=-bv(j-1);
FFF(j+1,j+1)=-bv(j-1);
end
FFF(62,63)=HV(61);
FFF(63,31)=-1;
for j=1:11,
FFF(62+j,63+j)=HN(j);
FFF(62+j,64+j)=-bn(j);
FFF(63+j,63+j)=-bn(j);
end
FFF(74,75)=HN(12);
FFF(75,30)=1;
FFF(75,32)=-1;
FFF(75,64)=-1;
xxx=FFF\F;
PE=xxx(1,1);
Unc=xxx(30,1);
Pnc=xxx(31,1);
for j=2:29,
gUv(j)=xxx(j,1);
end
for j=30:61,
gUv(j)=xxx(j+2,1);
end
for j=1:12,
gUn(j)=xxx(j+63,1);
end
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));
%/* refresh QLv, QCv, QLwv, u3v, QLn, QCn, QLwn, u3n,..., etc. */
gQLv(1) = 4.0*smpfrq*Ls(1)/2.0*leng(1)*gUv(1) - gQLv(1);
for j=2:NTS,
gQLv(j) = 4.0*smpfrq*(Ls(j-1)/2.0*leng(j-1) + Ls(j)/2.0*leng(j))*gUv(j) - gQLv(j);
end
gQLnts = 4.0*smpfrq*Ls(NTS)/2.0*leng(NTS)*Unc - gQLnts;
gQLv(NTS+1) = 4.0*smpfrq*Ls(NTS+1)/2.0*leng(NTS+1)*gUv(NTS+1) - gQLv(NTS+1);
for j=NTS+2:NING,
gQLv(j) = 4.0*smpfrq*(Ls(j-1)/2.0*leng(j-1) + Ls(j)/2.0*leng(j))*gUv(j) - gQLv(j);
end
gQLning = 4.0*smpfrq*Ls(NING)/2.0*leng(NING)*gUv(NING+1) - gQLning;
for j=1:NING,
gQCv(j)=4.0*smpfrq*Cp(j)*leng(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(0)+2) = 4.0*smpfrq*LNs(Nsin(0)+1)/2.0*NTslen*gUn(Nsin(0)+2) - gQLn(Nsin(0)+2);
end
if (NSIN == 2)
u3sin(1) = (gPsin(1)+gQLsin(1)-gVCsin(1))/bsin(1);
gQLsin(1) = 4.0*smpfrq*Lsin(0)*u3sin(1) - gQLsin(1);
gVCsin(1) = u3sin(1)/(smpfrq*Csin(0)) + gVCsin(1);
gQLns(1) = 4.0*smpfrq*LNs(Nsin(0))/2.0*NTslen*gUsin(1)-gQLns(1);
gQLn(Nsin(0)+2) = 4.0*smpfrq*LNs(Nsin(0)+1)/2.0*NTslen*gUn(Nsin(0)+2) - gQLn(Nsin(0)+2);
u3sin(2) = (gPsin(2)+gQLsin(2)-gVCsin(2))/bsin(2);
gQLsin(2) = 4.0*smpfrq*Lsin(1)*u3sin(2) - gQLsin(2);
gVCsin(2) = u3sin(2)/(smpfrq*Csin(1)) + gVCsin(2);
gQLns(2) = 4.0*smpfrq*LNs(Nsin(1))/2.0*NTslen*gUsin(2)-gQLns(2);
gQLn(Nsin(1)+2) = 4.0*smpfrq*LNs(Nsin(1)+1)/2.0*NTslen*gUn(Nsin(1)+2) - gQLn(Nsin(1)+2);
end
acoustic;
%/* 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=0:tnnum-1,
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=0:tnnum-1,
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 + -
显示快捷键?