cfbg_gaussianvisibility.m

利用传输矩阵法计算光纤啁啾brag光纤光栅滤波器的温度滤波特性

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  The following program is used to calculate propogation 
%  spectrum of a CFBG, given the original light source and
%  the expected spectrum, using transmiting matrix approach.
% 
% First  time:  2006-10-02
% Second time:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  

clear all; 
close all;
load sfs_all.mat;        % download the data of the original light source
L=0.08;                  % the length of a CFBG
T=300;                   % the initial temperature 
Neff=1.45;
lamda=[1528.08:0.16:1559.92]*1e-9;
M=200;                   % the number used in transmiting matrix appraoch
l=L/M;
C=200*(1e-9);            % linear chirped coefficient
%#####################################################  Gaussian function
tao=0.5;
G=sym('exp(-(t/1)^2)');
N1=20;
N2=190;
for i=1:N1
    g(i)=subs(G,-tao+(tao/(N1-1))*(i-1));
end
for i=(N1+1):200
    g(i)=1;
end
for i=(N1+1):N2
    g(i)=1;
end
for i=(N2+1):200
    g(i)=0.98*subs(G,(i-(N2+1))*tao/(200-(N2+1)));
end

%                micro-change   
for i=1:15
    g(i)=g(i)*0.9;
end
for i=45:55
    g(i)=g(i)*1.06;
end
for i=56:80
    g(i)=g(i)*1.09;
end
for i=110:200
    g(i)=g(i)*0.98;
end
%

v=0.41*g;


% obj=0.2*ones(1,200); % the expected spectrum is bandpass,with intensity of 0.2
% for i=1:6
%     obj(i)=0;
% end
% for i=163:200
%     obj(i)=0;
% end

% obj_x=1:200;
% obj_tao=50;
% obj_t=obj_x-100;
% obj=0.2*exp(-(obj_t/obj_tao).^2);

obj_x=1:200;
obj_tao=40;
obj_t=obj_x-40;
obj=1.16e-4*exp(-(obj_t/obj_tao).^2);


for i=1:M
%     r(i)=obj(i)/(sfs_all(i,1)/(0.1250e-3));
    r(i)=obj(i)/sfs_all(i,1);
    rr=r(i);
    lamda_B=lamda(i);
    f=sym('tanh(pi*x/((1+x/Neff)*lamda_B)*l)^2=rr'); % obtain delta_N from this equation       
    inter=solve(f);    % two solvtions are available here, and we only use the positive one later.
    b=subs(inter);
    delta_N(i)=b(1)/1.5;
    r_object(i)=tanh(pi*delta_N(i)/((1+delta_N(i)/Neff)*lamda_B)*l)^2;
end

 
count=0;
for i=1:M     % the transmitting approach                                                                                                                                                                                                                                                                                                                                                                                                                                                                                             
       F=[1,0;0,1];                                                                                                                                                                                                                                                                                                                                                                                                                                                                                               
    for k=1:M                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                         
	    Period(k)=lamda(k)/2/Neff;
        Beita(i,k)=2*pi*Neff/lamda(i);
        if k<2
           sigma=(Beita(i,k)-pi/Period(k))+2*pi*delta_N(k)/lamda(i); %(formula 17 from Erdogan's 'Fiber Grating Spectra' )
        else
           sigma=(Beita(i,k)-pi/Period(k))+2*pi*delta_N(k)/lamda(i)+8*pi*Neff^2*(l/2)/lamda(k)^2*(Period(k)-Period(k-1))/l;   
        end
        K=pi*v(k)*delta_N(k)/lamda(i);
        sub=sqrt(K^2-sigma^2);       
        F=F*[cosh(sub*L/M)-j*(sigma/sub)*sinh(sub*L/M),-j*(K/sub)*sinh(sub*L/M);j*(K/sub)*sinh(sub*L/M),cosh(sub*L/M)+j*(sigma/sub)*sinh(sub*L/M)];  %传输矩阵
        count=count+1
    end                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                     
        r_design(i)=(abs(F(3)/F(1)))^2;                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  
end
        initial_design=r_design;
        
for i=1:M
final(i)=(sfs_all(i,1)/(0.1250e-3))*r_design(i);
sfss(i)=sfs_all(i,1)/(0.1250e-3);
end

initial_final=final;
subplot(2,2,1);
plot(lamda,r_object)
title('object')
subplot(2,2,2);
plot(lamda,r_design)
title('design')
subplot(2,2,3);
plot(lamda,final)
title('final')
subplot(2,2,4);
plot(lamda,sfss)
title('sfs')

% plot(lamda,r_object);
% hold on;
% plot(lamda,r_design,'r');
% grid on;
% figure
% plot(lamda,final);
% grid on;

% subplot(2,1,1)
% plot(lamda,r_object);
% hold on;
% plot(lamda,r_design);
% grid on;
% subplot(2,1,2)
% plot(lamda,final);
% grid on;
% title('light source after filting by a CFBG');
% xlabel('Wavelength/m');        
% ylabel('Intensity');