soamodel.m

本人编写的半导体光放大器物理模型仿真程序。

% Simulate the SOA 
% Reference: JOS vol 21 No.9 Sept. 2004 Nielsen
%q=integrate(N,dz)
%dq=(q0-q)/Taosp-(G-1)Pin/hvA
%G=exp(Gama*a(q-NtrL))
%G0->q0
%L=500um, A=0.7*0.4um^2, Gama=0.6, I=70mA(20kA/cm^2)
%a1=5e-20m^2, Ntr=1e18cm^-3, Taosp=0.5ns,hv=0.945ev (1ev=1.602e-19J)
%Thegama: 1e16cm^-2, Gama*a1: cm^2, N*L*1e-4:cm^-2 
%Power: mw, P/(hvA):1/(16) 1e16ns^-1cm^-2
function SOADYNAMIC()
global Gama a1 Thegama0 Ntr Length Taosp G0 Esat hv Area t0 TimeResolution InputSignalSize
global Tnumber
% Parameter of SOA----------------------------
G0=16;%little signal gain dB
hv=0.945%ev, Wavelength=1.31;%um
a1=8.0;a2=0;a3=0;% a:1e-16cm^2, cm^-1 nm^-2,10^-17 cm^3nm
Length=300;%um
Area=0.4*0.4;%um^2
Gama=0.45;
Current=70;%mA
Ntr=2.6;%1e18cm^-3
Taosp=1;%ns

h=6.626;c=3;%10^-34J s , 10^8m/s
R=2.5;Tb=1/R;% rate 2.5Gb/s Tb:ns
TimeResolution=Tb/10;

InputCouplingLoss=0;OutputCouplingLoss=0;%Coupling loss

% Wave property-------------------------------
%sigal wave:rise cosin a=0.5; Gauss coeffiency:1
a=.5;Gauss=3;
ExtinctionRatio=10;%dB
DynamicRange=21;%dB
Attenuate2=29;%dB, Minimum Sensitiviy
Attenuate1=Attenuate2-DynamicRange;%dB
OnePower=0;%dBm, Logic 1
ZeroPower=OnePower-ExtinctionRatio;% Logic 0
%Input Time Domain Signal: dBm%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Signal=PBRS(200,31);% Pesudo Binary Randan Series
[time1,Pout1]=SignalGernerator(TimeResolution,Tb,a,OnePower,ZeroPower,Signal);
Signal=PBRS(200,73);
[time2,Pout2]=SignalGernerator(TimeResolution,Tb,a,OnePower,ZeroPower,Signal);
time=[time1 time2+TimeResolution+time1(size(time1,2))];
InputSignal=[Pout1-Attenuate1 Pout2-Attenuate2];%dBm
InputSignal=dBtA(InputSignal);%dBm to mw
t0=time(1);
InputSignalSize=size(InputSignal,2);

%Power Scan----------------------------------------------------
InputPower=-1*Attenuate2:0.5:-1*Attenuate1;%dBm
G0=dBtA(G0);
Thegama0=log(G0)/(Gama*a1)+Ntr*Length*1e-2;%cm^-2
%Thegama0=3+log(G0)/(Gama*a1);%cm^-2

%Esat=Psat*Tao0;%10^-12J
%Psat/Esat=1/ns
%-----------------------------------

InputPower=dBtA(InputPower-InputCouplingLoss);%dBm to mw
for i=1:size(InputPower,2)
%7.0    ThegamaEquibrium(i)=fzero(@(Thegama) ThegamaEquibriumfunc(Thegama,InputPower(i)),Thegama0);
    ThegamaEquibrium(i)=fzero(@ThegamaEquibriumfunc,Thegama0,[],InputPower(i));
end
GEquibrium=exp(Gama*a1*(ThegamaEquibrium-Ntr*Length*1e-2));

OutputPower=AtdB(InputPower.*GEquibrium);
InputPower=AtdB(InputPower);
% subplot(2,2,1),
%hold on;
plot(OutputPower,AtdB(GEquibrium),'b');%,InputPower,AtdB(GEquibrium));
%plot(InputPower,OutputPower);
xlabel('Output Power (dBm)');
ylabel('Gain (dB)');
grid on;
%---------------------------------------------------------------

% %Time Domain ---------------------------------------------------
% options = odeset('RelTol',1e-14,'AbsTol',[1e-14]);
% [time,ThegamaDynamic] = ode45(@ThegamaDynamicfunc,time,Thegama0,options,InputSignal);
% 
% GDynamic=exp(Gama*a1*(ThegamaDynamic-Ntr*Length*1e-2));
% 
% OutputSignal=AtdB(InputSignal'.*GDynamic);
% InputSignal=AtdB(InputSignal');
% 
% 
% %subplot(2,2,1),
% figure;
% plot(time,InputSignal);
% %subplot(3,1,2),plot(time,AtdB(GDynamic));
% xlabel('Time (ns)');
% ylabel('InputSignal (dBm)');
% 
% %subplot(2,2,2),
% figure;
% plot(time,OutputSignal);
% %subplot(3,1,2),plot(time,AtdB(GDynamic));
% xlabel('Time (ns)');
% ylabel('OutputSignal (dBm)');
% 
% SignalNum=length(OutputSignal);
% 
% %subplot(2,2,3),
% figure;
% EyeDiagram(OutputSignal(50:SignalNum/2),2);%
% 
% %subplot(2,2,4),
% figure;
% EyeDiagram(OutputSignal(SignalNum/2+50:SignalNum),2);%


%subfunction-------------------------------
% Equibrium function of Thegama
%dq=(q0-q)/Taosp-(G-1)Pin/hvA; dq=0

function f=ThegamaEquibriumfunc(Thegama,Pin)
global Gama a1 Thegama0 Ntr Length Taosp G0 Esat hv Area
GEquibrium=exp(Gama*a1*(Thegama-Ntr*Length*1e-2));
f=(Thegama-Thegama0)/Taosp+(GEquibrium-1)*Pin/(hv*Area*16);

% Dynamic function of Thegama
%dq=(q0-q)/Taosp-(G-1)Pin/hvA;

function f=ThegamaDynamicfunc(t,Thegama,Pin)
global Gama a1 Thegama0 Ntr Length Taosp G0 Esat hv Area t0 TimeResolution InputSignalSize
k=round((t-t0)/TimeResolution+1);
if k>InputSignalSize
   k=InputSignalSize;
end

GEquibrium=exp(Gama*a1*(Thegama-Ntr*Length*1e-2));
f=-(Thegama-Thegama0)/Taosp-(GEquibrium-1)*Pin(k)/(hv*Area*16);

function [y,time]=Pinc(xa)
%input signal wave
global Pc Tc Tb dt Signal t0 Snumber
time=t0:dt:-t0;
y=zeros(size(time));
for i=1:Snumber
   x=(time-(-round(Snumber/2)+i)*Tb)/Tc;
   y=y+exp(-x.^(2*xa)./2);
end
y=y*Pc;
%tadd=(ts0-Basen*dt):dt:(ts0-dt);
%time=[tadd,time];
%plot(time,y);

function y=Pin1(time)
%input control wave
global Ps
y=Ps.*ones(1,size(time,2));%super Gauss y=exp(-((t/T0).^(2*m))/2);

   
function y=dBtg(x)
global L
y=x./(1e-3.*L.*log10(exp(1)));

function y=gtdB(x)
global L
y=1e-3.*x.*L.*log10(exp(1));

function y=dBtA(x)
y=10.^(x/10);

function y=AtdB(x)
y=10.*log10(x);

function y=gtA(x)
global L
y=(exp(1)).^(x.*L*1e-4);

function y=g(ncarry,Lamda)
global Gama a1 n0 a2 Lam0 a3 ndc
y=Gama.*(a1.*(ncarry-n0).*100-a2.*1e6.*((Lamda-Lam0...
   +a3.*(ncarry-ndc).*0.01).^2));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Pesudo Randam Binary Sequence Generator (PRBS)
%r Generator Polynomial r Generator Polynomial 
% Sequence bits Generator Polynomial 
%7  [7 6 0] 26  [26 25 24 20 0] 
%8  [8 6 5 4 0] 27  [27 26 25 22 0] 
%10  [10 7 0] 29  [29 27 0] 
%11  [11 9 0] 30  [30 29 28 7 0] 
%12  [12 11 8 6 0] 31  [31 28 0] 
%15  [15 14 0] 34  [34 15 14 1 0] 
%16  [16 15 13 4 0] 35  [35 2 0] 
%Here we use 8bits z^8+z^6+z^5+z^4+1
%Seed:1~2^8-1

function y=PBRS(DigitalNumber,DigitalSeed)
y=[];
Feedback=[8 6 5 4 0]+1;
for i=1:DigitalNumber
    y=[y bitget(DigitalSeed,1)];
    FeedbackSum=mod(sum(bitget(DigitalSeed,Feedback)),2);
    DigitalSeed=bitshift(DigitalSeed,-1);
    DigitalSeed=bitset(DigitalSeed,Feedback(1),FeedbackSum);
end
    
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%55
%dt:     Time resolution
%Tsignal: Signal Period
%xa: Raise cosine coefficience
%Pmax Pmin
%DigitalSignal: Input pesudoramdon siganl
function [t,Pout]=SignalGernerator(dt,Tsignal,xa,Pmax,Pmin,DigitalSignal)
global Tnumber
x1=(-Tsignal/2:dt:0)/Tsignal;%1 up edge
x2=(0:dt:Tsignal/2)/Tsignal;% 1 down edge

%rise cosine func  xa:factor
Snumber=size(DigitalSignal,2);
Ps=1;
   BaseUpWave=sinc(x1).*cos(xa.*pi.*x1)./(1-4.*xa.*xa.*x1.*x1);
   BaseDownWave=sinc(x2).*cos(xa.*pi.*x2)./(1-4.*xa.*xa.*x2.*x2);
   PminLocal=min(BaseDownWave)*2;
   BaseUpWave=[-BaseDownWave+PminLocal BaseUpWave(2:size(BaseUpWave,2))];
   BaseEdgeNo=size(BaseUpWave,2)-1;
   Tnumber=BaseEdgeNo;
   for i=1:BaseEdgeNo
       BaseDownWave(i)=BaseUpWave(BaseEdgeNo+1-i);
   end
       BaseUpWave=BaseUpWave(2:BaseEdgeNo+1);
   
PminLocal=min(BaseDownWave);
  
%hold on;
   Pout=[];
   for i=2:Snumber-1
       if DigitalSignal(i)==1
           if DigitalSignal(i-1)==0
               Pout=[Pout BaseUpWave];
           else
               Pout=[Pout ones(1,BaseEdgeNo)];
           end
               

       else
           if DigitalSignal(i-1)==1
               Pout=[Pout BaseDownWave];
           else
               Pout=[Pout PminLocal*ones(1,BaseEdgeNo)];
           end
               
       end
   end
   PminLocal=min(Pout);
   Pout=Pmin+(Pmax-Pmin)/(1-PminLocal)*(Pout-PminLocal);
   t=linspace(0,dt*size(Pout,2),size(Pout,2));
   
function EyeDiagram(p,Tdisplay)
global Tnumber
hold off;
DataNumber=length(p);

PeriodNumber=round(DataNumber/Tdisplay/Tnumber);
Displaynumber=Tnumber*Tdisplay;

i=1;
plot(p(((i-1)*Displaynumber+1):i*Displaynumber));
hold on;
for i=2:PeriodNumber-1
    plot(p(((i-1)*Displaynumber+1):i*Displaynumber));
end
hold off;