on the synthesis of fiber bragg gratings by layer peeling.m

采用分离剥层算法以及傅立叶变化的方法

%%%%%参照文献《on the synthesis of fiber bragg grating by layer
%%%%%peeling》中第一个例子编写的程序,用自己编写的傅立叶程序变换
%%%%已试验成功，谨慎修改
clear all
clc    
Neff=1.446;
L=10;                       %%%光栅长10cm
N=500;                                       
delta=L/N;                   %%%分成500层
Lm=linspace(0,10,N);

M=1000;
step=4/(M-1);
wvl=[1548:step:1552]*1e-7;         %%%4nm的波长范围分成1000个点
wvlD=1550*1e-7;
deviate=2*pi*Neff*(1./wvl-1/wvlD); %%%失谐量
t=linspace(-20,20,M);
dw=157/(M-1);
dt=40/(M-1);

alpha=zeros(1,N);                  %%%用于存欲求的Rho
q=zeros(1,N);                      %%%用于存欲求的耦合系数                        
r=zeros(N,M);                      %%%用于存500个反射谱

r(1,:)=sqrt(0.9)*exp(-(deviate/19.2).^20);     %%%理想反射谱
%figure,plot(deviate,r(1,:))
%h=r*exp(i*deviate'*t)/(2*pi)*dw;              %%%傅立叶变换到时域
h1=r(1,:)*exp(i*deviate'*(t-10))/(2*pi)*dw;     %%%平移
win=sin(2*pi*t/40);                            %%%汉宁窗
h2=h1.*win;                             
h3=h2(M/2:M-1);                                %%%提取【0 10】的部分
h4=h3(N/4:3*N/4-1);
t3=linspace(0,20,M/2);
t4=linspace(5,15,N/2);
%figure,plot(t4,real(h4))
r(1,:)=h4*exp(i*t4'*deviate)*dt;                    %%%反变换到频域,求得物理可实现谱
%figure,plot(deviate,abs(r(1,:)))
%%%%%%求q
ph=exp(-i*2*(deviate)*delta);
for n=1:N
    rs=sum(r');                       %%%求和
    alpha(1,n)=rs(1,n)/M;             %%%文献中的(9)式              
    qabs=-atanh(abs(alpha(1,n)))/delta;             
    qphase=-phase(alpha(1,n));
    q(1,n)=qabs*exp(i*qphase);
    r(n+1,:)=ph.*((r(n,:)-alpha(1,n))./(1-conj(alpha(1,n))*r(n,:)));
end
figure,plot(Lm,real(q))                %%%%
%%%%%用传输矩阵法，求反射谱,参照文献《on the synthesis of fiber bragg grating by layer
%%%%%peeling》
for j=1:M
    Matrix_g=[1,0;0,1];
    for k=1:N
        kac=q(1,k);
        delta1=deviate(1,j);
        alpha=sqrt(abs(kac)^2-delta1^2);
        alphaL=alpha.*delta;

        T11=cosh(alphaL)+i*(delta1/alpha)*sinh(alphaL);
        T12=(kac/alpha)*sinh(alphaL);
        T21=(conj(kac)/alpha)*sinh(alphaL);
        T22=cosh(alphaL)-i*(delta1/alpha)*sinh(alphaL);
        Matrix_g=Matrix_g*[T11,T12;T21,T22];
    end
   TT11=Matrix_g(1);TT12=Matrix_g(3);TT21=Matrix_g(2);TT22=Matrix_g(4);
   reflect=-TT21/TT11;
   Reflectivity(j)=(abs(reflect)).^2;
   Q(j)=phase(reflect);
end
figure,plot(wvl,10*log10(Reflectivity))        %%%%画反射谱，纵轴单位为 db时，用这个方法比peeling中的好
%%%%%%%%求时延
P(1)=Q(1);P(2)=Q(2);P(3)=Q(3);    %%%求导就是求斜率
for l=4:M;
if(abs(Q(l-1)-Q(l))<=4.5)         %%%为什么是4.5？？？
    P(l)=((wvl(1))^2/(2*pi*3e8))*((Q(l-1)-Q(l))/(step*1e-9));   %%%%公式见《fiber grating spectra》（32）
else
    P(l)=((wvl(1))^2/(2*pi*3e8))*((Q(l-3)-Q(l-2))/(step*1e-9));
end
   temp_lamda(l-3)=wvl(l);
   temp_Delay(l-3)=-P(l);
end
Delay=10^12*temp_Delay;
figure,plot(temp_lamda,Delay)
%%%%%%用自己编写的程序，实现《on the synthesis of fiber bragg grating by layer
%%%%%%peeling》中第二个例子
%%%%已经试验成功，谨慎修改
clear all
clc  
Neff=1.446;
L=12;                       %%%光栅长10cm
N=600;                                       
delta=L/N;                   %%%分成500层
Lm=linspace(0,12,N);

M=1200;
step=4/(M-1);
wvl=[1548:step:1552]*1e-7;         %%%4nm的波长范围分成1000个点
wvlD=1550*1e-7;
deviate=2*pi*Neff*(1./wvl-1/wvlD); %%%失谐量
t=linspace(-24,24,M);
dw=157/(M-1);
dt=48/(M-1);

alpha=zeros(1,N);                  %%%用于存欲求的Rho
q=zeros(1,N);                      %%%用于存欲求的耦合系数                        
r=zeros(N,M);                      %%%用于存500个反射谱

r(1,:)=sqrt(0.9)*exp(-(deviate/12.3).^2).*exp(-i*0.13*deviate.^2);%%%理想反射谱,s=10.^12ps
figure,plot(wvl,abs(r(1,:)))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%未平移
%h=r(1,:)*exp(i*deviate'*(t))/(2*pi)*dw;                             %%%傅立叶变换到时域
%figure,plot(t,real(h))
%win_1=exp(-4*log(2)*(t).^2/(12));                                 %%%高斯窗
%win_1=0.5*(1+cos(2*pi*(t)/24));                                     %%%升余弦窗
%figure,plot(t,win_1)
%h_1=h.*win_1;
%figure,plot(t,real(h_1))
%r(1,:)=h_1*exp(i*(t)'*deviate)*dt;
%figure,plot(deviate,abs(r(1,:)))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%平移
h1=r(1,:)*exp(i*deviate'*(t-12))/(2*pi)*dw; 
%figure,plot(t,real(h1))
win=exp(-4*log(2)*(t-12).^2/(8));                   %%%汉宁窗
%win=0.5*(1+cos(2*pi*(t-12)/12));                   %%%%升余弦窗
h2=h1.*win;                             
h3=h2(M/2:M-1);                                      %%%提取【0 10】的部分
h4=h3(N/4:3*N/4-1);
t3=linspace(0,24,M/2);
%figure,plot(t3,real(h3))
%t4=linspace(6,18,N/2);
%figure,plot(t4,real(h4))
r(1,:)=h4*exp(i*t4'*deviate)*dt;                     %%%反变换到频域,求得物理可实现谱
figure,plot(wvl,abs(r(1,:)))
%%%%%求q
ph=exp(-i*2*(deviate)*delta);
for n=1:N
    rs=sum(r');                       %%%求和
    alpha(1,n)=rs(1,n)/M;             %%%文献中的(9)式              
    qabs=-atanh(abs(alpha(1,n)))/delta;             
    qphase=-phase(alpha(1,n));
    q(1,n)=qabs*exp(i*qphase);
    r(n+1,:)=ph.*((r(n,:)-alpha(1,n))./(1-conj(alpha(1,n))*r(n,:)));
end
figure,plot(Lm,abs(q))                %%%%
%win_q=exp(-4*log(2)*(Lm-6).^2/(8));
%figure,plot(Lm,win_q)
%q1=q.*win_q;
%figure,plot(Lm,abs(q1))
%%%%%%%求反射谱
for j=1:M
    Matrix_g=[1,0;0,1];
    for k=1:N
        kac=q(1,k);
        delta1=deviate(1,j);
        alpha=sqrt(abs(kac)^2-delta1^2);
        alphaL=alpha.*delta;

        T11=cosh(alphaL)+i*(delta1/alpha)*sinh(alphaL);
        T12=(kac/alpha)*sinh(alphaL);
        T21=(conj(kac)/alpha)*sinh(alphaL);
        T22=cosh(alphaL)-i*(delta1/alpha)*sinh(alphaL);
        Matrix_g=Matrix_g*[T11,T12;T21,T22];
    end
   TT11=Matrix_g(1);TT12=Matrix_g(3);TT21=Matrix_g(2);TT22=Matrix_g(4);
   reflect=-TT21/TT11;
   Reflectivity(j)=(abs(reflect)).^2;
   Q(j)=phase(reflect);
end
figure,plot(deviate,10*log10(Reflectivity))        %%%%画反射谱，纵轴单位为 db时，用这个方法比peeling中的好%%
%%%%%%%%%%%%%%%求时延
P(1)=Q(1);P(2)=Q(2);P(3)=Q(3);    %%%求导就是求斜率
for l=4:M;
if(abs(Q(l-1)-Q(l))<=4)         %%%为什么是4.5？？？
    P(l)=((wvl(1))^2/(2*pi*3e8))*((Q(l-1)-Q(l))/(step*1e-9));   %%%%公式见《fiber grating spectra》（32）
else
    P(l)=((wvl(1))^2/(2*pi*3e8))*((Q(l-3)-Q(l-2))/(step*1e-9));
end
   temp_lamda(l-3)=wvl(l);
   temp_Delay(l-3)=-P(l);
end
Delay=10^8*temp_Delay;
figure,plot(temp_lamda,Delay)
axis([1.5495e-4 1.5505e-4 0 1000])