dg.m

本微波光子发生器仿真软件实现了对一种基于两个光纤光栅的新光子微波发生器的仿真分析

clear;
clc;
%一一一一一一一一一一一一一参数说明
%光源部分
w0=1550.1;        %THz，即中心波长为 154On。
dw=3;%标准高斯函数中的一个参数，单位THz，
dw_3dB=dw*2*log(2);%%光源3dB带宽
dw=dw_3dB/(2*log(2));
Sta_w=w0-10;
End_w=w0+10;
Step_w=0.01;
%w=Sta_w:step_w:End_w;%考虑的频率范围(注意是角频率)
%MZI部分
%dww=1;   %MZI的透过率周期为1THz
w=Sta_w:Step_w:End_w;         %THz
S=exp(-((w-w0)./dw).^2);%光源的高斯包络
ss=plot(w,S);
set(ss,'Color',[.7,.1,.9]);
xlabel('波长（nm）');
ylabel('幅度');
title('光源的高斯包络');
figure;
tp=1;
%switch tp
 %   case 1
y1=[1540:0.01:1560];
%yd=1550; 
neff=1.48;
s=1;
dneff=0.0002;
 l=1*1e6;                         %光栅长度
yd1=1550;                        %FBG1中心波长
k1=pi*s*dneff./y1;
s1=2*pi*neff*(1./y1-1/yd1);
 d1=2*pi*dneff./y1;
 ds1=s1+2*pi*dneff./y1;
 w1=(k1*l).^2-(ds1*l).^2;
r1=6*(sinh(w1.^0.5)).^2./((cosh(w1.^0.5)).^2-ds1.^2./k1.^2);%FBG1反射率
yc1=(1+dneff/neff)*yd1;
FWHM1=s*dneff*sqrt(1+(yd1/s*dneff*l)^2)/(4*neff);   %FBG1半波带宽
sd1=plot(y1,r1);
%sd=plot(y,r);
set(sd1,'Color',[.7,.1,.9]);
xlabel('波长(nm)');
ylabel('反射率');
title('FBG1的反射谱');
figure;
r1=r1.*S;
ss1=plot(y1,r1);
set(ss1,'Color',[.7,.1,.9]);
xlabel('反射波长(nm)');
ylabel('功率(mW)');
title('FBG1取得的反射波');
figure;
%hold on;
%均匀FBG1反射谱的取得==================================================
neff2=1.68;
y2=[1540:0.01:1560];
yd2=1549.000;
dneff=0.0002;
neff=1.48;
s=1;
k2=pi*dneff./y2;
l=1*1e6;   
s2=2*pi*neff*(1./y2-1/yd2);
d2=2*pi*dneff./y2;
ds2=s2+2*pi*dneff./y2;
w2=(k2*l).^2-(ds2*l).^2;
r2=6*(sinh(w2.^0.5)).^2./((cosh(w2.^0.5)).^2-ds2.^2./k2.^2);
yc2=(1+dneff/neff)*yd2;
FWHM2=s*dneff*sqrt(1+(yd2/s*dneff*l)^2)/(4*neff);
r2=r2.*S;
sd2=plot(y2,r2);
%sd1=plot(y1,r1);
%sd=plot(y,r);
set(sd2,'Color',[.1,.1,.9]);
xlabel('波长(nm)');
ylabel('反射率');
title('FBG2的反射谱');
figure;
r2=r2.*S;
ss2=plot(y1,r2);
set(ss2,'Color',[.7,.1,.9]);
xlabel('反射波长(nm)');
ylabel('功率(mW)');
title('FBG2取得的反射波');
figure;    
coupled=r1+r2;%耦合波
y=[1540:0.01:1560];
sd3=plot(y,coupled);
%sd1=plot(y1,r1);
%sd=plot(y,r);
set(sd3,'Color',[.1,.1,.3]);
xlabel('波长(nm)');
ylabel('功率(mW)');
title('FBG1和FBG2反射波的耦合波');
figure;
%=================以上为均匀光纤光栅反射谱取得和耦合部分=====================
%case 2
%啁啾光纤光栅的反射谱的取得。
%clear;
%clc;
C=0.17e-9;       %啁啾
L=140;           %长度
n0=0.0;          %直流调制
n11=0.5;          %交流调制 
n=1.4568;        %有效折射率
period0=534.5;   %周期
span=0.3;        %谱范围
p=0.0003;
c=2.99792458e8;
L=L*(1e-3);
n0=n0*(1e-4);
n11=n11*(1e-4);
period0=period0*(1e-9);
span=span*(2e-9);
p=p*(1e-9);
M=140;
dz=L/M;
centerwavelength=2*(n+n0)*period0;
wavelength=(centerwavelength-span/2):p:(centerwavelength+span/2);
M1=length(wavelength);
hh=waitbar(0,'计算中...');
for j=1:1:M1
  waitbar(j/M1,hh)  
v=wavelength(j);
period1=period0-C/2;
n1=n11*exp((-(4*log(2)*(-L/2)^2))/((L/3)^2));  % for Gaussian apodization---q=0
k=pi*n1/v;
h=2*(n+n0)*pi/v-pi/period1;
omig=sqrt(k*k-h*h);
T(1,1)=cosh(omig*dz)-(i*h*sinh(omig*dz))/omig;
T(1,2)=-i*k*sinh(omig*dz)/omig;
T(2,1)=i*k*sinh(omig*dz)/omig;
T(2,2)=conj(T(1,1));
for q=1:1:M
    B=T;
    period(q)=period0+(+q/M-0.5)*C;
     n1(q)=n11*exp((-(4*log(2)*(q*L/M-L/2)^2))/((L/3)^2));  % for Gaussian apodization
    k=pi*n1(q)/v;
    h(q)=2*(n+n0)*pi/v-pi/period(q);
    omig(q)=sqrt(k*k-h(q)*h(q));
    T(1,1)=cosh(omig(q)*dz)-(i*h(q)*sinh(omig(q)*dz))/omig(q);
    T(1,2)=-i*k*sinh(omig(q)*dz)/omig(q);
    T(2,1)=i*k*sinh(omig(q)*dz)/omig(q);
    T(2,2)=conj(T(1,1));
    T=T*B;
%以下增加了一个随机的相差
 %^^^^^^^^^^^^^^^^^^^
  anf0=rand(1)*10;
    anf=fix(anf0)*pi/180/30;
     phase(1,1)=exp(-i*anf/2);
    phase(1,2)=0;
    phase(2,1)=0;
    phase(2,2)=exp(i*anf/2);
T=phase*T;
% ^^^^^^^^^^^^^^^^^^^^^   
end
    A=T*[1;0];
    r(j)=A(2,1)/A(1,1);
    transmission(j)=1/abs(A(1,1));
    transmission1(j)=abs(A(1,1));
    transmission2(j)=abs(A(2,1));
    Reflectivity(j)=(abs(r(j)))^2;
    sita(j)=angle(r(j));
end    
%figure;
%plot(wavelength*(1e9),sita, 'r');
figure;
ss3=plot(wavelength*(1e9),Reflectivity, 'r');
set(ss3,'Color',[.7,.1,.9]);
xlabel('波长 nm');
ylabel('反射率')
title('啁啾光纤光栅反射谱');
axis auto;
diudao1=diff(sita)/p;
diudao1(M1)=diudao1(M1-1);
for j=1:1:M1
    v=wavelength(q);
    delay(j)=-v^2*diudao1(j)/(2*pi*c);
end
if abs(delay(1))>1e-9
    delay(1)=0;
end
for j=2:1:M1
    if delay(j)>7e-9
        delay(j)=delay(j-1);
    end
    if delay(j)<-7e-12 
        delay(j)=delay(j-1);
    end
end
%figure;
%plot(wavelength,delay*(1e12));
%title('时延');
axis auto;
xlabel('wavelength ');
ylabel('Delay (ps)');
for j=1:1:M1
    a(j,1)=wavelength(j);
    a(j,2)=delay(j);
end
P1=0; P2=0; P3=0; P4=0;
for i=1:1:M1
    P1=P1+a(i,1);
    P2=P2+a(i,2);
    P3=P3+a(i,1)^2;
    P4=P4+a(i,1)*a(i,2);
end
avedispersion=(P4*M1-P1*P2)/((P3*M1-P1^2))
b=(P3*P2-P1*P4)/(M1*P3-P1*P1)
for j=1:1:M1
    dispersion(j)=avedispersion;
    ripple(j)=delay(j)-(avedispersion*wavelength(j)+b);
end
figure;
ss4=plot(wavelength,dispersion); 
set(ss4,'Color',[.7,.1,.9]);
title('色散');
axis auto;
xlabel('波长（nm） ');
ylabel('色散 x1000ps/nm');
%figure;
%plot(wavelength,ripple*(1e12));
%title('纹波');
axis auto;
xlabel('wavelength ');
ylabel('时延纹波 (ps)');
%figure;
%plot(wavelength,Reflectivity,wavelength,delay*(1e9),wavelength,ripple*(1e9)+1);
hold off;
figure;
plotyy(wavelength,Reflectivity,wavelength,delay*(1e9))
figure;
ss5=plot(wavelength,transmission);
set(ss5,'Color',[.7,.1,.9]);
title('透射谱');
close(hh)
%otherwise
 %       error('error!tp=1->NRZ;2->RZ;3->Manchester');
%end
%=======================================================-==================
%参数定义--=======================================================--外差探测部分
%光电检测器仿真子程序
%[amp,x]=photo_detect(pt)
%pt->input array
%unit:nA
nt=size(coupled,2);
lambda_c=1.51*10^-6;
%dt=10^-16;
dt=1.25e-7;
h=6.626*10^-34;                   %普朗克常量
c=3*10^8;                         %光速常量
e=1.602*10^-19;                   %库伦常数
vs=c/lambda_c;                    %光频率
yita=0.8;                         %光电转换效率
Id=5e-9;
figure;
for ii=1:length(r1)
    aa1(ii,1)=y1(ii);
    aa1(ii,2)=r1(ii);
end
for jj=1:length(r1)
    if aa1(jj,2)==max(r1)
        marky=jj;
    end
    if aa1(jj,1)==1549
        markx1=jj;
    end
    if aa1(jj,1)==1551
        markx2=jj;
    end
  end
%yy=[1549:0.001:1551];
for ki=markx1:markx2
    %r1(ki)=r1(ki+1)-r1(ki);
    r1(ki)=r1(ki+5)-r1(ki-12);
end
alfa=e*yita/(h*vs);
iout=e*yita*r1/(h*vs);
f=[1:0.0145:30];
sat1=plot(f,iout);
set(sat1,'Color',[.7,.1,.9]);
xlabel('频率（GHz）');
ylabel('幅度（A.U）');
title('FBG1反射波经过PD检测输出的光电流');
%=======================================================-==================
figure;
for iii=1:length(r2)
    aa2(iii,1)=y2(iii);
    aa2(iii,2)=r2(iii);
end

for jjj=1:length(r2)
    if aa2(jjj,2)==max(r2)
        marky2=jjj;
    end
    if aa2(jjj,1)==1549
        markx12=jjj;
    end
    if aa2(jjj,1)==1551
        markx22=jjj;
    end
  end
%yy=[1549:0.001:1551];
for kij=markx12:markx22
    %r1(ki)=r1(ki+1)-r1(ki);
    r2(kij)=r2(kij+5)-r2(kij-17);
end
alfa=e*yita/(h*vs);
iout2=e*yita*r2/(h*vs);
f=[1:0.0145:30];
sat3=plot(f,iout2);
set(sat3,'Color',[.7,.1,.9]);
xlabel('频率（GHz）');
ylabel('幅度（A.U）');
title('FBG2反射波经过PD检测输出的光电流');
figure;
ioutL=iout+iout2;
ss6=plot(f,ioutL);
set(ss6,'Color',[.7,.1,.9]);
xlabel('频率（GHz）');
ylabel('幅度（A.U）');
title('耦合波经过PD检测输出的低频光电流部分');
figure;
for iiii=1:length(coupled)
    aa2(iiii,1)=y2(iiii);
    aa2(iiii,2)=coupled(iiii);
end

for jjjj=1:length(coupled)
    if aa2(jjjj,2)==max(coupled)
        marky22=jjjj;
    end
    if aa2(jjjj,2)==max(coupled)
        marky21=jjjj;
    end
    if aa2(jjjj,1)==1547
        markx12=jjjj;
    end
    if aa2(jjjj,1)==1552
        markx22=jjjj;
    end
  
end
%yy=[1549:0.001:1551];
for kiij=markx12:markx22
    %r1(ki)=r1(ki+1)-r1(ki);
    coupled(kiij)=coupled(kiij+5)-coupled(kiij-21);
end
alfa=e*yita/(h*vs);
iout3=e*yita*coupled/(h*vs);
f0=linspace(60,80,2001);
ss7=plot(f0,iout3*10);
%sat2=plot(f,iout3);
set(ss7,'Color',[.7,.1,.9]);
xlabel('频率（GHz）');
ylabel('幅度（A.U）');
title('耦合光经过PD检测输出的高频光电流部分');
figure;
f3=linspace(15,80,4002);
for mm=500:length(coupled)
    ioutL(mm-499)=ioutL(mm);
end
for mn=1:2*length(coupled)
    if mn<=2001
        IOUT(mn)=ioutL(mn);
    end
    if mn>2001
        IOUT(mn)=iout3(mn-2001)*2.5;
    end
end
ss8=plot(f3,IOUT);
set(ss8,'Color',[.7,.1,.9]);
%sat2=plot(f,iout3);
xlabel('频率（GHz）');
ylabel('幅度（A.U）');
title('耦合光经过PD检测输出的光电流');
figure;
%高频带通滤波器部分
Order=1;
CenterFrequency=65;
%中心频率CenterFrequency单位为GHz
Bandwidth=0.01;%单位GHz
Factor=((f3-CenterFrequency)/(Bandwidth/2.0)).^(2.0*Order); 
TransferFunction=exp(-0.5*0.693147180559945309417*Factor);
%滤波器传输函数
RF=TransferFunction.*IOUT;
for i=1:length(f3)
    if RF(i)<0
        RF(i)=-RF(i);
    end
end
ss9=plot(f3,RF);
set(ss9,'Color',[.7,.1,.9]);
xlabel('频率（GHz）');
ylabel('幅度（A.U）');
title('经高频带通滤波器输出的毫米波信号');
%高频带通滤波器部分结束
%总程序结束=======================================================-==================