shixi.m

晶体光滤波器的设计。用matlab完成

%=======定义变量x========%
syms x;
%======A5,B5以及系数=====%
MA(1,6)=-0.0310-0.0695*x+0.0972*x^2+0.3174*x^3-0.5658*x^4+0.2525*x^5;       %A5
ma(6,:)=[-0.0310 -0.0695 0.0972 0.3174 -0.5658 0.2525];                     %A5的系数
MA(2,6)=-0.0310*x^5-0.0695*x^4+0.0972*x^3+0.3174*x^2-0.5658*x+0.2525;       %A_5
A55=x^5-MA(1,6)*MA(2,6);  
%==========求根==========%
root=solve(A55);                                                          
%======求得root的根======% 
root =[
[-.96777138118745328374426985895375+.42624722068466751003525435674718*i]
[-.96777138118745328374426985895375-.42624722068466751003525435674718*i]
[-.86541979253568298369398365624900+.38116727614038152624718392667623*i]
[-.86541979253568298369398365624900-.38116727614038152624718392667623*i]
[ .30442290514249513302749028847461+.16654664561981403800124920640494*i]
[ .30442290514249513302749028847461-.16654664561981403800124920640494*i]
[ 2.5281965662491176594826252516403+1.3831503821278832658354182950700*i]
[ 2.5281965662491176594826252516403-1.3831503821278832658354182950700*i]
[                                                                    -1]
[                                                                    -1]];
%======选择其中五个根======% 
MB(1,6)=(x-root(3))*(x-root(4))*(x-root(5))*(x-root(6))*(x-root(9));        %B5
expand(MB(1,6));
mb(6,:)=[9.2871    19.7071   10.0563  -3.4846   -2.1209    1 ];                        %B5的系数
%========求修正系数=======% 
suma=0;sumb=0;
for t=1:6;
    suma=suma+ma(6,t)^2;                                                    %A5的系数平方相加
    sumb=sumb+mb(6,t)^2;                                                    %B5的系数平方相加
end 
alpha=sqrt((1-suma)/sumb);
MB(1,6)=alpha*MB(1,6);                                                      %多项式修正
mb(6,:)=alpha*mb(6,:);                                                      %系数修正   
expand(MB(1,6));                                                            %修正后多项式展开
seita(6)=-atan(mb(6,6)/ma(6,6));
%=====循环求A4,A3,A2,A1,B4,B3,B2,B1以及系数======% 
for m=5:-1:1;j=m;
    MA(1,m)=x*( cos(seita(j+1))*MA(1,m+1)-sin(seita(j+1))*MB(1,m+1) );
    expand(MA(1,m));
    ma(m,:)=cos(seita(j+1))*ma(m+1,:)-sin(seita(j+1))*mb(m+1,:);
    for n=1:5;
        ma(m,n)=ma(m,n+1);
    end 
        ma(m,m+1)=0;
    MB(1,m)=sin(seita(j+1))*MA(1,m+1)+cos(seita(j+1))*MB(1,m+1); 
    expand(MB(1,m));
    mb(m,:)=sin(seita(j+1))*ma(m+1,:)+cos(seita(j+1))*mb(m+1,:);
    seita(j)=-atan(mb(m,j)/ma(m,j));                                        %求角度    
end
%=====循环求相对角θ======%     
fai(1)=0;
for k=2:7;
    fai(k)=seita(k-1)+fai(k-1);
end
%=======验证琼斯矩阵===========%
J=1;   
for s=2:6;
    J=J*[cos(fai(s)) -sin(fai(s));sin(fai(s)) cos(fai(s))]*[x 0; 0 1]*[cos(fai(s)) sin(fai(s));-sin(fai(s)) cos(fai(s))];
end 
s=7;
J=J*[cos(fai(s)) -sin(fai(s));sin(fai(s)) cos(fai(s))];
expand(J);
%=====φ,θ换算成角度====%
    fai=180/pi*fai;
    seita=180/pi*seita;
%===========绘图============%     
   c=3e8;n=0.2;f=5e10;L=c/(n*f); 
   wavemin=1548.7*10^(-9);
   wavemax=1549.3*10^(-9);
   wavelength=wavemin:0.001*10^(-9):wavemax;
   deltafai=n*L*2*pi./wavelength;
%===========理想情况下============%   
   y=-0.0310-0.0695*exp(-i*deltafai)+0.0972*exp(-2*i*deltafai)+0.3174*exp(-3*i*deltafai)-0.5658*exp(-4*i*deltafai)+0.2525*exp(-i*5*deltafai);
   z=20*log10(abs(y));
%===========实际情况下============%  
   yy=-0.0916-0.0738*exp(-i*deltafai)+ 0.1403*exp(-2*i*deltafai)+0.3100*exp(-3*i*deltafai)-0.5794*exp(-4*i*deltafai)+0.2547*exp(-i*5*deltafai);
   zz=20*log10(abs(yy)); 
   plot(wavelength,z,'g',wavelength,zz,'r');   
   axis([wavemin,wavemax,-55,1]);