ctr.m

基于半导体激光器的matlab程序产生混沌

%chaos transmitter and receiver 
%without Langevin noise 

function ctr 
    clc
    clear
    syms t w;  
T1=clock;
  %激光器参数
    p.c=3*10^8;p.h=6.63*10^-34;
    p.V=1.5*10^-16;            %active region volume
    p.q=1.6*10^-19;            %电荷所带电量
    p.Gn=2.125*10^-12;         %增益斜率
    p.af=5.5;                  %线宽增强因子
    p.gm=1*10^-5;            %自发辐射因子
    p.eps=3*10^-23;         %增益饱和参量
    p.yt=0.4;                  %光场限制因子
    p.tn=2*10^-9;           %载流子寿命
    
    p.tl=9*10^-12;
    p.n0=4*10^23;              %透明载流子密度
    p.nt=9.9*10^23;            %阈值载流子密度
    p.th_e=0;
    
    p.tpt=2*10^-12;          %光子寿命
    p.tpr=2*10^-12;          %光子寿命
   
    %*****控制项--电流--********************
    p.f=2*10^9;
            ith=0.012;
            rb=1.8;
            md=0.0;
    p.it=rb*ith;
    p.ip=md*abs(rb-1)*ith;
    
    p.ir=p.it;
    
    %*****控制项--反馈--********************
            R2=0.3;
            p.wt=2*pi*p.c/(1.55*10^-6);  p.wr=2*pi*(p.c/(1.55*10^-6)+0*10^9);
            Ret=1.0/100;                 Rer=1.0/100;
            let=7.50/100;                ler=7.50/100;
    p.kft=(1-R2)*sqrt(Ret/R2);           p.kfr=(1-R2)*sqrt(Rer/R2);
    p.tet=2*let/p.c;                     p.ter=2*ler/p.c;
    
    Rj=30/100;
    p.kj=(1-R2)*sqrt(Rj/R2); 
    p.tj=1*10^-9;
            
            kj=p.kj
            kt=p.kft
            kr=p.kfr
            lagt=p.tet;                        %反馈光注入延迟时间
            lagr=p.ter
    %*****dde23************************           
            fs=4*10^10;              
            ts=0; 
            tf=2*10^-8;
            xt0=[9.9*10^23;4.02*10^20;0];         %t<=0时的值
            xr0=[9.9*10^23;4.02*10^20;0];
%_____________________________________________________________
    
            solt=dde23(@sctrt,lagt,xt0,[ts,tf],[],p);
            
            solr=dde23(@sctrr,lagr,xr0,[ts,tf],[],p,solt);
         %________________________________________________________________    

                    t=1*tf/5:1/fs:tf;
                    %***********************************
                    [yt,ypt]=deval(solt,t-p.tj);
                    Pt=p.V*p.h*(p.wt/(2*pi)).*yt(2,:)/(2*p.yt*p.tpt);       %光功率
                    [ft,pwt]=spectrum(fs,Pt);

                    %***********************************
                    [yr,ypr]=deval(solr,t);
                    Pr=p.V*p.h*(p.wr/(2*pi)).*yr(2,:)/(2*p.yt*p.tpr);       %光功率
                    [fr,pwr]=spectrum(fs,Pr);
                            
           figure(1)
                    plot(ft,pwt,'black')
                    hold on
                    plot(fr,pwr,'r')
                    xlabel('frequency (Hz)');ylabel('power dBm');
                    title('spectrum')
           figure(2)
                    plot(t*10^9,Pt,'black','LineWidth',1.3)
                    hold on
                    plot(t*10^9,Pr,'r','LineWidth',1.3)
                    xlabel('time t/ns');ylabel('output power p/W');
           figure(3)
                    subplot(211),plot(yt(1,:)/p.nt,Pt,'black')
                    subplot(212),plot(yr(1,:)/p.nt,Pr,'black')
                    xlabel('n/nth');ylabel('output power P/W');
           figure(4)
                    plot(Pt,Pr,'*black','MarkerSize',2)
                    xlabel('output power Pt/W');ylabel('output power Pr/W');
                    %eyediagram(P,6,1*10^-10)
           figure(5)
                    plot(t*10^9,yt(3,:),'black','LineWidth',1.5)
                    hold on
                    plot(t*10^9,yr(3,:),'r','LineWidth',1.5)
                    xlabel('time t/ns');ylabel('phase');
           figure(6)
                    plot(yt(3,:),yr(3,:),'*black','MarkerSize',2)
                    xlabel('phase of transmitter');ylabel('phase of receiver');
                    %eyediagram(P,6,1*10^-10)       
                    %*********
                    std(Pr-Pt)
                    dp=mean(abs(Pt/max(Pt)-Pr/(max(Pr))))
                    pht=abs(yt(3,:));phr=abs(yr(3,:));
                    dph=mean(abs(pht/max(pht)-phr/(max(phr))))
                    mean((Pr-Pt-mean(Pr-Pt)).^2)
                    %*********
                    t=2*tf/5:1/fs:tf*0.8;
                    y1=deval(solr,t);
                    y2=deval(solr,t-1*tf/5);
                    figure(7)
                    [c,lag]=xcorr(y1(2,:),y2(2,:),'coeff');
                    plot(10^9*lag/fs,c)
                    
    fprintf('\n');
   % fprintf('\t');
    fprintf('CPUTIME IS');
    lll=etime(clock,T1);
    disp(lll);
%**************************************************************************
%**************************************************************************