soa_halfwave00.m

用MATLAB编程实现半导体光放大器（SOA）的ASE噪声特性分析

%SOA速率方程
clear;
clc;
y=0.892;      % Molar fraction of Arsenide in the active region
L=0.7e-3;     %有源层长度 m 
d=0.4e-6;     %有源层高度 m
W=0.4e-6;       %有源层宽度 m
T=0.45;       %限制因子
Kg=0.9e-10;   % bandgap shrinkage coefficient
n1=3.22;      %有源区反射率
c1=2.998e8;   %真空中光速 m/s
Arad=1e7;          %线性复合系数 s^-1 &
Brad=5.6e-16;      %双分子复合系数 m^3/s
Anrad=3.5e8;
% Bnrad=0e-16;
Caug=3.0e-41;      %耦合系数   m^6/s
a=1.35;         %Bandgap energy quadratic coefficient
b=-0.775;       %Bandgap energy quadratic coefficient
c=0.149;        %Bandgap energy quadratic coefficient
me=4.10e-32;    %effective mass of electron in the CB
mhh=4.19e-31;   %effective mass of a heavy hole in the VB
mhl=5.06e-32;   %effective mass of a light hole in the VB
R1=5e-5;        %放大器输入端面的功率反射系数
R2=5e-5;        %放大器输出端面的功率反射系数
r1=sqrt(R1);
r2=sqrt(R2);
neq0=3.22;
dneq=-1.34e-26;
K0=6200;
K1=7500e-24;
h=6.626e-34;    %普朗科常数 J*s
Tem=300;        %绝对温度
k0=1.38e-23;     %玻耳兹曼常数
e=1.602e-19;    %电子电量 c
vsig=c1/1.55e-6;   %信号频率
p0=1e24;        %割线法的初始点
p1=4e24;
number=0;       %循环次数
ng=3.55;             %sqrt((3.22^2-3.167^2)*(T/(2-T))+3.167^2)
vv=c1/(2*ng*L); %自发辐射谱的频率间隔   2e11
p=10;
M=40;           %放大器分为10段
N=35;           %自发辐射谱分为10段
l=L/M;          %每小段长度
etaout=0.5;
etain=0.5;


    %单独计算第M+1段的载流子      
    Psigin=10^(-30/10-3);                 %输入探测功率 W
    Ssigp=ones(1,M+1)*0;
    Ssign1=ones(1,M+1)*0;
    Ssign2=ones(1,M+1)*0;
    Esign=ones(1,M+1)*0;
    Esigp=ones(1,M+1)*0;
           
    Sspn=ones(M+1,N)*0;     %建立正向和负向传播的 (M+1)*N 阶数组来存储自发辐射谱
    Sspp1=ones(M+1,N)*0;
    Sspp2=ones(M+1,N)*0;
    
    for s=1:1:N
        v(s)=1.8738e+014+(s-1)*vv*p;         %自发辐射频率1.2308e-019/h   载流子为1e24时的频率
        lambda(s)=c1/v(s);
    end
    
NUMtol=90;        

%对容忍度进行判断，决定是否进行新的循环
while NUMtol>0 
    for s=1:1:N
        Sspp1(1,s)=(Sspp1(1,s)+Sspp2(1,s))/2;
        Sspp2(1,s)=Sspp1(1,s);         
    end
    
    Ssign1(1)=(Ssign1(1)+Ssign2(1))/2;
    Ssign2(1)=Ssign1(1);
        
    Ssigp(1)=(Psigin*etain)/(h*vsig)*(1-R1)+R1*Ssign1(1);
    Esigp(1)=sqrt(Ssigp(1));
     
    n(1)=secant1(p0,p1,Ssigp(1),Ssign1(1),Sspp1(1,:),Sspn(1,:));    %割线法求解载流子密度
   
    for i=1:1:M %计算第M段到第1段的载流子密度   
%         neq(i)=neq0+dneq*n(i);
        beta(i)=1.4512e+007;%2*pi*neq(i)*vsig/c1;    %计算传输系数
        ain(i)=K0+T*K1*n(i);            %计算每段的内部损耗
        Esigp(i+1)=Esigp(i)*exp((0.5*(T*gm(n(i),vsig)-ain(i)))*l)*exp((-beta(i)*1i)*l);
        Ssigp(i+1)=abs(Esigp(i+1))^2;
        for s=1:1:N  
          Rsp(i,s)=T*gm(n(i),v(s))*vv*p*fc(n(i),v(s))*(1-fv(n(i),v(s)))/(fc(n(i),v(s))-fv(n(i),v(s)));  %2表示两个正交方向
          Sspp1(i+1,s)=(Sspp1(i,s)+Rsp(i,s)/(T*gm(n(i),v(s))-ain(i)))*exp((T*gm(n(i),v(s))-ain(i))*l)-Rsp(i,s)/(T*gm(n(i),v(s))-ain(i));
%           Sspp1(i+1,s)=Sspp1(i,s)*exp((T*gm(n(i),v(s))-ain(i))*l)+Rsp(i,s)/(T*gm(n(i),v(s))-ain(i))*(exp((T*gm(n(i),v(s))-ain(i))*l)-1);
        end
        n(i+1)=secant1(p0,p1,Ssigp(i+1),Ssign1(i+1),Sspp1(i+1,:),Sspn(i+1,:));
    end
 
    Ssigout=Ssigp(M+1)*(1-R2);
    Esign(M+1)=r1*Esigp(M+1);
    Ssign1(M+1)=abs(Esign(M+1))^2;
    Sspn(M+1,:)=R1.*Sspp1(M+1,:);  
    Psigout=Ssigout*h*vsig*etain;
    %前半个循环结束
    
    %后半个循环开始
    for j=M:-1:1
        Esign(j)=Esign(j+1)*exp((0.5*(T*gm(n(j),vsig)-ain(j)))*l)*exp((-beta(j)*1i)*l);
        Ssign1(j)=abs(Esign(j))^2;
        for s=1:1:N
          Sspn(j,s)=(Sspn(j+1,s)+Rsp(j,s)/(T*gm(n(j),v(s))-ain(j)))*exp((T*gm(n(j),v(s))-ain(j))*l)-Rsp(j,s)/(T*gm(n(j),v(s))-ain(j));
        end
    end
     
        
    tolsig=abs((Ssign1(1)-Ssign2(1))/compare(Ssign1(1),Ssign2(1)));  %判断反向信号是否达到容忍度
    if  tolsig>1e-3
         NUMsig=1; 
    else NUMsig=0;
    end    
        
    for s=1:1:N                   %判断正向噪声是否达到容忍度
        Sspp1(1,s)=Sspn(1,s)*R2;
        tolsp(s)=abs((Sspp1(1,s)-Sspp2(1,s))/compare(Sspp1(1,s),Sspp2(1,s)));
    end     
    NUMsp=0;                      %记录有几个容忍度大于1e-4
    for s=1:1:N
         if tolsp(s)>1e-3
            NUMsp=NUMsp+1; 
         end
    end      
    NUMtol=NUMsp+NUMsig
    number=number+1
    tolsp(N)
    if number==50
        break,
    end
 end

 
 %画图
 G1=zeros(1,N);
 G2=zeros(1,N);
 for s=1:1:N
     vares(s)=v(s)-vv/2;
     for i=1:1:M
         G1(s)=G1(s)+(T*gm(n(i),v(s))-ain(i))*l;
         G2(s)=G2(s)+(T*gm(n(i),vares(s))-ain(i))*l;
     end
     Gres(s)=(1-R1)*(1-R2)*exp(G1(s))/(1-sqrt(R1*R2)*exp(G1(s)))^2;
     Gares(s)=(1-R1)*(1-R2)*exp(G2(s))/(1+sqrt(R1*R2)*exp(G2(s)))^2;
     r(s)=(4*sqrt(R1*R2)*exp(G1(s)))/(1-sqrt(R1*R2)*exp(G1(s)))^2;
     K(s)=1/sqrt(1+r(s)^2);
     sigmaN(s)=2*etaout*K(s)*Sspp1(M+1,s)*(1-R2)/(vv*p);          %噪声光子速率谱密度
     sigmaASE(s)=sigmaN(s)*h*v(s);             %噪声功率谱密度*Gares(s)/Gres(s)
     Pase(s)=sigmaASE(s)*1e-10*v(s)^2/c1;                   %噪声功率
     PdBm(s)=10*(log10(Pase(s))+3);      
 end
 
Gsig=10*log10(Psigout/Psigin);
NF=10*log10(sigmaASE(11)/(h*vsig*Psigout/Psigin)+etaout/(Psigout/Psigin));
%NF=10*log10(sigmaASE(11)/(h*vsig*Gsig)+etaout/Gsig);  
%NF=10*log10(sigmaASE(11)/(h*vsig*Psigout/Psigin)+etaout/(Psigout/Psigin))+1;    %噪声数dB  

plot(lambda*1e9,Pase)
xlabel('Wavelength (nm)'); 
ylabel('pase (w)');

figure
plot(lambda*1e9,sigmaASE)
xlabel('Wavelength (nm)'); 
ylabel('sigmaASE(w)');

figure
plot(lambda*1e9,PdBm) 

xlabel('Wavelength (nm)'); 
ylabel('Power (dBm)');



for s=1:1:N
   G(s)=10*log10(Sspp1(M+1,s)/Sspp1(1,s));  
   r(s)=(4*sqrt(R1*R2)*exp(G(s)))/(1-sqrt(R1*R2)*exp(G(s)))^2;
   K(s)=1/sqrt(1+r(s)^2);
end   

Ssppall=0;
Sspnall=0;
for i=1:1:N
    Ssppall=Ssppall+2*Sspp1(M+1,i);
    Sspnall=Sspnall+2*Sspn(1,i);
end

figure
for i=1:1:M+1
    carrier(i)=n(i);
    length(i)=l*(i-1)*1e6;
end
plot(length,carrier)
xlabel(' Length (\mum)'); 
ylabel('Carrier Density (m^-^3)');


figure
for i=1:1:N  
    gg(i)=gm(n(20),v(i));
end
plot(lambda*1e9,gg)
xlabel('Wavelength (nm)'); 
ylabel('gg(m^-^1)');

figure
plot(lambda*1e9,NF)      %输入信号VS噪声系数
xlabel('Wavelength (nm）'); 
ylabel('Noise Figure (dB)');
 for i=1:1:15
  
     NF1(i)=NF(i);
 end
 plot(lambda*1e9,NF1)
 xlabel('Wavelength (nm）'); 
 ylabel('Noise Figure (dB)');


figure
subplot(211)
plot(length,Ssppall(1,:)/1e15)
hold on
plot(length,Sspnall(1,:)/1e15)
xlabel('SOA Length (\mum)'); 
ylabel('ASE photon rate (10^1^6s^-^1)');

subplot(212)
plot(length,Ssppall(13,:)/1e15)
hold on
plot(length,Sspnall(13,:)/1e15)
xlabel('SOA Length (\mum)'); 
ylabel('ASE photon rate (10^1^5s^-^1)');
 


 

% G1=ones(1,N)*0;
% G2=ones(1,N)*0;
% for s=1:1:N
%     vres(s)=1.2308e-019/h+(s-1/2)*vv;         %自发辐射频率1.2308e-019   载流子为1e24时的频率
%     vares(s)=vres(s)-vv/2;
%     for i=1:1:M
%         G1(s)=G1(s)+(T*gm(n(i),vres(s))-ain(i))*l;
%         G2(s)=G2(s)+(T*gm(n(i),vares(s))-ain(i))*l;
%     end
%     Gres(s)=(1-R1)*(1-R2)*exp(G1(s))/(1-sqrt(R1*R2)*exp(G1(s)))^2;
%     Gares(s)=(1-R1)*(1-R2)*exp(G2(s))/(1+sqrt(R1*R2)*exp(G2(s)))^2;
%     NF(s)=10*log10(sigmaASE(s)/(h*v(s)*Gres(s))+0.5/Gres(s));
% end
% 
% figure
% plot(lambda*1e9,NF)