extra_pole1.rar_coupling matrix_mxc_极点_耦合滤波器_耦合矩阵

clear all; N=4; RL=22; tz=[1.3217,1.8082]; pw=poly(tz); ppw=conv(pw,pw); lz=length(tz); for k1=(lz+1):N tz(k1)=Inf; end %求U,V多项式 syms w; U1=w-1/tz(1); V1=((w^2-1)^0.5)*((1-1/(tz(1)^2))^0.5); Ulast=U1; Vlast=V1; for k2=1:(N-1) U=NextU(Ulast,Vlast,tz(k2+1)); V=NextV(Ulast,Vlast,tz(k2+1)); Ulast=U; Vlast=V; end; %求解E(S),P(S),F(S)的表达式 uw=sym2poly(Ulast); fz=roots(uw); fw=poly(fz); ffw=conv(fw,fw); rip=(1/(10^(RL/10)-1))^0.5*abs(polyval(pw,1)/polyval(fw,1)); eew=[zeros(1,length(ffw)-length(ppw)),ppw]+ffw*(rip^2); tzee=roots(eew); tze=tzee(find(imag(tzee)>0)); pwr=pw/rip; ew=poly(tze); ps=poly(j*tz); es=poly(j*tze); fs=poly(j*fz); if mod(N-lz,2)==0 ps=j*ps; end % 求[ABCD]矩阵多项式 efs=es+fs; %只需使用AS和BS,若要CS和DS，则用ES-FS efs0=es-fs; efs1=zeros(1,N+1); %初始化 efs10=zeros(1,N+1); efs2=zeros(1,N+1); efs20=zeros(1,N+1); for k1=(N+1):-2:1 % 偶次项 0，2，4 efs1(k1)=j*(imag(efs(k1))); efs10(k1)=j*(imag(efs0(k1))); efs2(k1)=real(efs(k1)); efs20(k1)=real(efs0(k1)); end for k2=N:-2:1 %奇次项 efs1(k2)=real(efs(k2)); efs10(k2)=real(efs0(k2)); efs2(k2)=j*(imag(efs(k2))); efs20(k2)=j*(imag(efs0(k2))); end if lz==N epr = rip/sqrt(rip^2-1); msl = epr/rip/(epr+1); else epr = 1.0; msl = 0.0; end y21n=ps/rip; %对AS，BS赋值 if mod(N,2)==0 % N为偶,BS为复奇 as=efs1; ds=efs10; bs=efs2; cs=efs20; else % N为奇，BS为复偶 as=efs2; ds=efs20; bs=efs1; cs=efs10; end if lz==N y21n=y21n-j*msl*bs; end y21n=j*y21n; % ps/rip %%%%%%% 提取极点综合 %第一步：移除第一个极点，提取相位长度theta1;并求解[A'B'C'D']; s01=j*tz(2); theta1=atan(polyval(bs,s01)/(j*polyval(ds,s01))); as1=as*cos(theta1)-j*cs*sin(theta1); bs1=bs*cos(theta1)-j*ds*sin(theta1); cs1=cs*cos(theta1)-j*as*sin(theta1); ds1=ds*cos(theta1)-j*bs*sin(theta1); %as1,bs1除以（s-s01)得到中间多项式axs和bxs ss0=[1,-s01]; axs=deconv(as1,ss0); bxs=deconv(bs1,ss0); laxs=length(axs); lbxs=length(bxs); if laxs==N+1 else for k1=1:N+1-laxs axs=[0,axs]; end for k1=1:N+1-lbxs bxs=[0,bxs]; end end %计算并联谐振器对的留数b01; b01=polyval(ds1,s01)/polyval(bxs,s01); %计算得到c''(s)和d''(s),p''(s) cs20=cs1-b01*axs; ds20=ds1-b01*bxs; as2=axs; bs2=bxs; cs2=deconv(cs20,ss0); ds2=deconv(ds20,ss0); y2n=deconv(y21n,ss0); lcs2=length(cs2); lds2=length(ds2); if lcs2==N+1 else for k1=1:N+1-lcs2 cs2=[0,cs2]; end for k1=1:N+1-lds2 ds2=[0,ds2]; end end %%提取第二个极点 s02=j*1.3217; theta12=atan(polyval(bs2,s02)/(j*polyval(ds2,s02))); as1=as2*cos(theta12)-j*cs2*sin(theta12); bs1=bs2*cos(theta12)-j*ds2*sin(theta12); cs1=cs2*cos(theta12)-j*as2*sin(theta12); ds1=ds2*cos(theta12)-j*bs2*sin(theta12); %as1,bs1除以（s-s02)得到中间多项式axs和bxs ss0=[1,-s02]; axs=deconv(as1,ss0); bxs=deconv(bs1,ss0); laxs=length(axs); lbxs=length(bxs); if laxs==N+1 else for k1=1:N+1-laxs axs=[0,axs]; end for k1=1:N+1-lbxs bxs=[0,bxs]; end end %计算并联谐振器对的留数b01; b01=polyval(ds1,s02)/polyval(bxs,s02); %计算第二次得到c''(s)和d''(s),p''(s) cs20=cs1-b01*axs; ds20=ds1-b01*bxs; as2=axs; bs2=bxs; cs2=deconv(cs20,ss0); ds2=deconv(ds20,ss0); y2n=deconv(y2n,ss0); %%第三次提取 s03=j*1e100; theta23=atan(polyval(bs2,s03)/(j*polyval(ds2,s03))); %Inf,Matlab无法计算，用1e300 lcs2=length(cs2); lds2=length(ds2); if lcs2==N+1 else for k1=1:N+1-lcs2 cs2=[0,cs2]; end for k1=1:N+1-lds2 ds2=[0,ds2]; end end as3=as2*cos(theta23)-j*cs2*sin(theta23); bs3=bs2*cos(theta23)-j*ds2*sin(theta23); cs3=cs2*cos(theta23)-j*as2*sin(theta23); ds3=ds2*cos(theta23)-j*bs2*sin(theta23);