Continue to Site

Welcome to EDAboard.com

Welcome to our site! EDAboard.com is an international Electronics Discussion Forum focused on EDA software, circuits, schematics, books, theory, papers, asic, pld, 8051, DSP, Network, RF, Analog Design, PCB, Service Manuals... and a whole lot more! To participate you need to register. Registration is free. Click here to register now.

Frequency Comb generation in Matlab

Status
Not open for further replies.

Sur76

Newbie level 1
Joined
May 26, 2011
Messages
1
Helped
0
Reputation
0
Reaction score
0
Trophy points
1,281
Activity points
1,359
Hello,

I am designing a optical comb to be used in DWDM systems, can anyone please help me in giving some conceptual inputs regarding how the laser can be gain switched to achieve the comb. I have a code that solves the rate equations to give photon density, carrier density and phase of the light emmitted. So can you please throw some light on this topic as I'm a bit new to this domain.

Thanking you in advance.
Sur76



P.S: The following is the code I have.

clear all


ind2 = 1;
steps = 0:0.1:20;

% Bind = 1;
% Bias_step = 16:4:24;
%
% for Biasing = Bias_step
for ind = steps

%% %%%%%%%%%%%%%%%%%%%%% Setup %%%%%%%%%%%%%%%%%%%%%%%%%%%%
Nbits = 1024*10;
freq = ind*1e9;
%freq = 0.1*1e9;
w = 2*pi*freq;
% bitrate = 9e10;
bitrate = 300e9;
dt = 1/bitrate;
df = 1/dt;
time = (-Nbits/2:Nbits/2-1)*dt;
FREQ=(-(Nbits/2):(Nbits/2-1))*1/(dt*Nbits);
FREQGHz = FREQ.*1e-9;

%% %%%%%%%%%%%%%%%% Drive Current %%%%%%%%%%%%%%%%%%%%%%%%%%

% amp = 10e-3; %mA
%
% Bias = 50e-3; %mA

%%%%default%%%%
amp = 4e-3; %mA
Bias = 50e-3; %mA
%%%%%%%%%%%%%%%
I = (amp/2)*sin(2*pi*freq*time);
%Bias = Biasing*1e-3;
I = I + Bias;

%I = Bias.*ones(size(time));

% I = I - mean(I);



%% %%%%%%%%%%%%%%%%%%%% Parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%
g0 = 1e-12; % Differential gain coeff. !!!!1.6e-12
Nom = 1.4e23; % Transparency density (m^-3) 1.4e23
V = 11e-17; % Volume of active layer 11e-17
tp = 2e-11; % Photon lifetime (s) !!!!2.4e-11
tn = .3e-9; % Carrier recombination lifetime !!!!.03e-8
OCF = 0.35; % Mode confinement factor 0.35
B = 0.0000; % Beta, spontaneous emission factor 0.0
q = 1.6e-19; % Charge of electron (C) 1.6e-19
alpha = 6.8; % Linewidth Enhancement Factor
%Ac = 8e-3; % Amplitude of carriers
Ac = amp;
%I_bias = 70e-3; % Bias Current
I_bias = Bias;
deltaf = -01e9; % Frequency Detuning -11
deltaw = 2*pi*deltaf; %Convert frequency to radians
Si = 0e20; %Injection Level (Injected photon density) - change to 0.7e20
Kc = 2.5e11; %Injected light coulping coefficient
if Si == 0;
deltaw = 0;%it there's no injection then ignore the detuning in the phase equation
end
h = 6.625e-34; % Plancks constant
R = 0.32; % Reflectivity in cavity
c = 3e8; % Speed of light
n = 3.63; % Refractive index
Ar = .03e-12; % Area of the active region
lamda = 1550e-9; % Wavelength of output light
f_laser = c/lamda;% Frequency of the output light
responsivity = 0.6; % Responsivity of the detector
fm1 = 1e7; %Data Modulation Rate
bitperiod = 1/fm1; %Data bit period Period = 1/f
numcarriers = 5; %number of subcarriers
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% THIS SECTION GETS THE STEADY STATE VALUES FOR PHASE PHOTON DENSITY AND CARRIER % % % % DENSITY IT THEN USES THEM TO WORK OUT THE SMALL SIGNAL RESPONSE OF THE LASER WHICH % % IT PLOTS THIS IS THE MODULATION RESPONSE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% the coefficients for a quartic equation in N0 to determine the steady
% state carrier density
% a_var*x^4 + b_var*x^3 + c_var*x^2 + d_var*x + e_var = 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
a_var = OCF^2*V^2*q^2*g0^2*tp^2*(alpha^2 + B^2 - 2*B + 1);
b_var = OCF*V*q*g0*tp*(I_bias*OCF*g0*tn*tp*(alpha^2 - B + 1) ...
+ V*q*(Nom*OCF*g0*tp*(alpha^2 + (B-1)^2) + alpha^2 ...
+ alpha*2*tp*deltaw -B + 1));
c_var = (I_bias^2*OCF^2*g0^2*tn^2*tp^2*(alpha^2 + 1) ...
+ 2*I_bias*OCF*V*q*g0*tn*tp*(2*Nom*OCF*g0*tp*(alpha^2 - B + 1) + 2*alpha^2 ...
+ 4*alpha*tp*deltaw - B + 2) + V^2*q^2*(4*Kc^2*Si*g0*tn*tp^2 ...
+ Nom^2*OCF^2*g0^2*tp^2*(alpha^2 + B^2 - 2*B + 1) ...
+ 2*Nom*OCF*g0*tp*(alpha^2 + 2*alpha*tp*deltaw - B + 1) + alpha^2 ...
+ 4*alpha*tp*deltaw + 4*deltaw^2*tp^2 + 1));
d_var = tn*(I_bias^2*OCF*g0*tn*tp*(Nom*OCF*g0*tp*(alpha^2 + 1) + alpha^2 ...
+ 2*alpha*tp*deltaw + 1) + I_bias*V*q*(2*Kc^2*Si*g0*tn*tp^2 ...
+ Nom^2*OCF^2*g0^2*tp^2*(alpha^2 - B + 1) ...
+ Nom*OCF*g0*tp*(2*alpha^2 + 4*alpha*tp*deltaw - B + 2) ...
+ alpha^2 + 4*alpha*tp*deltaw + 4*deltaw^2*tp^2 + 1) ...
+ 2*Kc^2*Nom*Si*V^2*q^2*g0*tp^2);
e_var = I_bias^2*tn^2*(Nom^2*OCF^2*g0^2*tp^2*(alpha^2 + 1) ...
+ 2*Nom*OCF*g0*tp*(alpha^2 + 2*alpha*tp*deltaw +1) + alpha^2 ...
+ 4*alpha*tp*deltaw + 4*deltaw^2*tp^2 + 1) ...
+ 4*I_bias*Kc^2*Nom*Si*V*q*g0*tn^2*tp^2;
%Get the roots of the above equation
N0_eq = [a_var -2*b_var c_var -2*d_var e_var];
N0_Roots = roots(N0_eq);
N0 = N0_Roots(4); %N0, the steady state carrier density 4!!
% S0 the steady state photon density is worked out by rearranging the
% carrier density rate equation
S0 = ((I_bias*tn)- (N0*q*V))/(tn*g0*(N0-Nom)*q*V);
% Injection Ratio
ratio = Si/S0;
% Output power is worked out using Photon Density
Pout = S0.*c./(2*OCF*n)*h*f_laser*Ar*(1-R);
%To work out steady state phase, Phi0, both equations give same answer but different %sign
if Si == 0;
Phi0 = 0;
else
% Phi0 = acos(((S0/tp) - (OCF*g0*(N0-Nom)*S0))/(2*Kc*sqrt(S0*Si)));
Phi0 = asin((-alpha/(2*tp) - (deltaw) + ...
((alpha/2)*OCF*g0*(N0-Nom)))/(Kc*sqrt(Si/S0)));
end
% Set up frequency points to plot the frequency response
% freqM = [10e6:10e6:20e9];
% wm = 2*pi*freqM;
%
% % %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % % Using the small signal derivation of the rate equations to get the modulation %response Modulation response is given as change in photon number with respect to % %change in input current Input current is Ac and change in photon density is s1
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% X = 2*Kc*sqrt(Si*S0)*cos(Phi0);
% Y = Kc*sqrt(Si/S0)*sin(Phi0);
% O = (j.*wm + (X/(2*S0)));
% P = (j.*wm + (1/tn) + g0*S0);
% if Si == 0
% s1 = -((Ac/(q*V))*g0*S0*OCF)./(wm.^2 - (wm).*(i*g0*S0) + (wm).*(i/(tn)) - (g0*S0)/tp);
% end
% if Si > 0
% s1 = (((OCF*Ac*g0*S0)./(P.*q*V)).*(1 - Y*alpha./O)) ... %top line
% ./((j.*wm) + (X/(2*S0)) + (Y^2./O) + (X*Y*alpha*g0./(P.*O)) - ...
% (S0*Y*alpha*g0./(P.*O*tp)) - (X*g0./P) + (g0*S0./(P.*tp)));
% end
%% %%%%%%% Plot Resonance Freq. %%%%%%%%%
% figure(2);
% s = abs(s1); %Get absolute values of the change in photon density
% logs = 20*log10(s/Ac); %Get the log of change in photon number with respect to change
% %in input current (The modulation response)
% norms = logs - logs(1); %Normalise it to the first value
% plot(freq.*1e-9,norms); % and plot
% axis([0 20 -30 40]);
% title('Modulation Response');
% xlabel('Frequency (GHz)'); ylabel('Response (dB)');
% hold on; grid;


%% %%%%%%%%%%%%%%%%%%% Initial Parameters %%%%%%%%%%%%%%%%%%%%%
S0 = abs(S0);
N0 = abs(N0);
Phi0 = abs(Phi0);

P1 = S0;
N1 = N0;
Phi1 = Phi0;
% P1 = 1e23;
% N1 = Nom;
% Phi1 = -3;

%% %%%%%%%%%%%%%%%%%%%%%% Calculate Rate Equations %%%%%%%%%%%%%%%%%
dt=time(2)-time(1);

N=N1*ones(size(time));
P=P1*ones(size(time));
Phi=Phi1*ones(size(time));

for k=2:length(time),
dPdt= ((OCF*g0*(N(k-1)-Nom) - (1/tp))*P(k-1)) + ((OCF*B*N(k-1))/tn)...
+ ((2*Kc)*sqrt(Si*P(k-1))*cos(Phi(k-1)));
dNdt= (I(k-1)/(q*V)) - (g0*(N(k-1)-Nom)*P(k-1)) - (N(k-1)/tn);


P(k)=dPdt*dt+P(k-1); %Photon Density
N(k)=dNdt*dt+N(k-1); %Carrier Density

if P(k) <= 0
P(k) = 100;
end

if Si == 0
Phi(k) = 0;
else
dPhidt = (alpha/2)*(OCF*g0*(N(k-1)-Nom) - 1/tp) - (deltaw) - (Kc*sqrt(Si/P(k-1))*sin(Phi(k-1)));
Phi(k)=dPhidt*dt+Phi(k-1);
end

end

%% %%%%%%%%%%%%%%%%%%%% Get Power %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Power_leftside = P.*c/(2*OCF*n)*h*f_laser*Ar*(1-R);
Power_leftside = Power_leftside*1e3; %mW
%Power_leftside = Power_leftside - mean(Power_leftside);
Power_leftsidef=abs(fftshift(fft(fftshift(Power_leftside))));

%Note: Taking FFT includes transcients at start => high power
%Normal for frequency response graph to start at 0dBm for small
%signal analysis


% Sig_Pwr = max(Power_leftside(1,(length(Power_leftside))/2:length(Power_leftside)));
% response(ind2) = Sig_Pwr; %mW

Sig_Pwr = max(Power_leftside(1,(length(Power_leftside)-5000):length(Power_leftside)))...
- min(Power_leftside(1,(length(Power_leftside)-5000):length(Power_leftside))); % Power signal amplitude after transcients mW
response(ind2) = 10*log10(Sig_Pwr); %dBm
%response(ind2) = 10*log10(Sig_Pwr./(amp.^2))
Sig_Count(ind2) = Sig_Pwr;
ind2 = ind2 + 1;
end
% Mult_Res(Bind,:) = response;
% Bind = Bind + 1;
% ind2 = 1;
% end

%% %%%%%%%%%%%%%%%%%%%%%% Graphs %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% figure(3)
% plot(FREQGHz, 10*log10(Power_leftsidef/max(Power_leftsidef)))
% title('Power v Frequency')
% xlabel('Frequency (GHz))')
% ylabel('Power(dB)')
% axis([min(FREQGHz) max(FREQGHz) -60 0])

response_graph = steps;

figure, plot(response_graph, response)
title('Laser Frequency response')
xlabel('Frequency (GHz)')
ylabel('Power (dBm)')
hold on; grid;

% figure(6)
% plot(10*log10(abs(fftshift(fft(fftshift(P(length(P)/2:length(P))))))))
%
%
% ax=[0 max(FREQGHz) -60 0];
% If=abs(fftshift(fft(fftshift(I)))).^2;


% figure(1)
% plot (time,I.*1e3)
% title('Tx time signal')
% xlabel('time')
% ylabel('mA')
%
%
% figure(2)
% plot(FREQGHz,10*log10(If/max(If)))
% title('Tx spectrum')
% axis(ax)
 
Last edited:

Status
Not open for further replies.

Part and Inventory Search

Welcome to EDABoard.com

Sponsor

Back
Top