From 55db794162103c66c165456dd1cb18e6687ee3f4 Mon Sep 17 00:00:00 2001 From: Silas Oettinghaus Date: Mon, 15 May 2023 10:32:24 +0200 Subject: [PATCH] WIP --- Classes/Amplifier.m | 2 +- Classes/Fiber.m | 125 ++++++++++++++++++++++++++++++++++++++++--- Classes/Photodiode.m | 63 ++++++++++++++++++++++ comm_tb.m | 16 +++--- 4 files changed, 190 insertions(+), 16 deletions(-) create mode 100644 Classes/Photodiode.m diff --git a/Classes/Amplifier.m b/Classes/Amplifier.m index d11eaee..dc90618 100644 --- a/Classes/Amplifier.m +++ b/Classes/Amplifier.m @@ -97,7 +97,7 @@ classdef Amplifier if obj.amp_mode == "output_power" seemsright = pow_out_dbm == obj.amplification_db; - else obj.amp_mode == "gain" + elseif obj.amp_mode == "gain" seemsright = pow_out_dbm == pow_in_dbm+ obj.amplification_db; end diff --git a/Classes/Fiber.m b/Classes/Fiber.m index aca6995..30ff317 100644 --- a/Classes/Fiber.m +++ b/Classes/Fiber.m @@ -3,6 +3,7 @@ classdef Fiber % Detailed explanation goes here properties + fsimu fiber_length alpha D @@ -10,6 +11,11 @@ classdef Fiber lambda0 gamma dphimax + + b2 + b3 + alpha_lin + linstep end methods @@ -17,6 +23,7 @@ classdef Fiber %FIBER Construct an instance of this class % Detailed explanation goes here arguments + options.fsimu options.fiber_length = 0 options.alpha = 0.2 options.D = 17 @@ -26,13 +33,21 @@ classdef Fiber options.dphimax = 5e-3 end - obj.fiber_length = options.fiber_length*1000; %km - obj.alpha = options.alpha; - obj.D =options.D*1e-6; - obj.Dslope =options.Dslope; - obj.lambda0 =options.lambda0; - obj.gamma =options.gamma; - obj.dphimax =options.dphimax; + obj.fsimu = options.fsimu; + obj.fiber_length = options.fiber_length*1000; %km + obj.alpha = options.alpha; + obj.D = options.D*1e-6; + obj.Dslope = options.Dslope*1e3; + obj.lambda0 = options.lambda0*1e-9; + obj.gamma = options.gamma; + obj.dphimax = options.dphimax; + + obj.b2 = -obj.D*obj.lambda0^2/(2*pi*Constant.LightSpeed); + obj.b3 = ((obj.lambda0.^2/(2*pi*Constant.LightSpeed)).^2*obj.Dslope); + obj.alpha_lin = obj.alpha/10*log(10)/1000; + + + end @@ -40,12 +55,108 @@ classdef Fiber %METHOD1 Summary of this method goes here % Detailed explanation goes here + N = length(opt_in); + faxis = linspace(-obj.fsimu/2,obj.fsimu/2,N+1); + faxis = ifftshift(faxis(:,1:end-1)); + faxis = faxis'; + + obj.linstep = -obj.alpha_lin/2 - 2*1j*pi^2*obj.b2*faxis.^2 - 4/3*1j*pi^3*obj.b3*faxis.^3; + + + opt_out = obj.NLSE(opt_in); + %attenuate nase + end + + function [yout] = NLSE(obj, xin) + + maxPow = obj.gamma.*max(abs(xin).^2); + + Leff = obj.dphimax / maxPow ; + + dz = Leff; + + z_prop = 0; + + yout = fft(xin); + + yout = ((yout).*exp(obj.linstep*dz/2)); + + while true + + yout = ifft(yout); + + Leff = dz; + + power = abs(yout).^2; + Hnl = exp( -1j*obj.gamma*power*Leff); + yout = yout .* Hnl; + + z_prop = z_prop + dz; + + maxPow = obj.gamma*max(abs(yout).^2); + Leff = obj.dphimax/maxPow; + + dz_new = Leff; + + if z_prop + dz_new > obj.fiber_length + dz_new = obj.fiber_length - z_prop; + break + end + + yout = fft(yout); + + yout = ((yout).*exp(obj.linstep*(dz/2+dz_new/2))); + + dz = dz_new; + + end + + yout = fft(yout); + + yout = ((yout).*exp(obj.linstep*(dz/2+dz_new/2))); + + yout = ifft(yout); + + Leff = dz; + + power = abs(yout).^2; + + Hnl = exp( -1j*obj.gamma*power*Leff); + + yout = yout .* Hnl; + + yout = fft(yout); + + yout = ((yout).*exp(obj.linstep*(dz/2))); + + yout = ifft(yout); + end + + function yout = process_nonlinstep(xin,stepsize) + + end + + function yout = process_linstep(xin,stepsize) end end end + + + + + + + + + + + + + + diff --git a/Classes/Photodiode.m b/Classes/Photodiode.m new file mode 100644 index 0000000..1378061 --- /dev/null +++ b/Classes/Photodiode.m @@ -0,0 +1,63 @@ +classdef Photodiode + %PGOTODIODE Summary of this class goes here + % Detailed explanation goes here + + properties + fsimu + responsivity + dark_current + temperature + + end + + methods + function obj = Photodiode(options) + %PHOTODIODE Construct an instance of this class + % Detailed explanation goes here + arguments + options.fsimu + options.responsivity = 1; + options.dark_current = 0; + options.temperature = 20; + end + + obj.fsimu = options.fsimu; + obj.responsivity = options.responsivity; + obj.dark_current = options.dark_current; + obj.temperature = options.temperature; + + end + + function yout = process(obj,xin) + %METHOD1 Summary of this method goes here + % Detailed explanation goes here + + k = Constant.Boltzmann; + T = obj.temperature + 273.15 ; %celsius + 273 = kelvin + R = 50; %resistance of phdiode (50ohm is typical value) + + % Magnitude squared detection + yout = sum( abs(xin) .^2*obj.responsivity,1) ; + + % Shot Noise + shot_noise = sqrt(k * obj.fsimu .* yout) .* randn(1,size(yout,1)); + + yout = yout + shot_noise; + + % Thermal Noise + therm_current_psd = (2 * k * T / R ) ; %squared + + Bw = obj.fsimu; %is this correct? shouldnt it be the bandwidth of the actual component? e.g. 70GHz? + therm_noise_pow = therm_current_psd * 2 * Bw; %squared + + therm_noise = sqrt(therm_noise_pow) .* randn(1,size(yout,1)); + + yout = yout + therm_noise; + + % Dark Current + yout = yout + obj.dark_current ; + + end + end +end + diff --git a/comm_tb.m b/comm_tb.m index 785399e..dd876b6 100644 --- a/comm_tb.m +++ b/comm_tb.m @@ -16,7 +16,6 @@ kover = 16; fsimu = kover * fdac ; - %SIMULATE for i = 1:log2(M) @@ -28,7 +27,6 @@ pamData = pam_mapping(bitpattern,M,0); shapedData = applyPulseShaping(pamData,fsym,fdac); - awg = AWG('preset','M8196A','fdac',fdac,'kover',kover,'lpf_active',1); awgSignal = awg.process_channel(shapedData); @@ -43,7 +41,12 @@ laserfield = eml.process(filtered); att = Amplifier("amplification_db",10,"amp_mode","gain","type","ideal","saturation_mode",0,'saturation_power',10); att_out = att.process(laserfield); -%fib = Fiber("fiber_length",2,"alpha",0.2,"D",17,"lambda0",1550); +fib = Fiber("fsimu",fdac*kover,"fiber_length",20,"alpha",0.2,"D",17,"lambda0",1550); +fib_out = fib.process(att_out); + +phdiode = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20); +phdiod_out = phdiode.process(fib_out); + figure; @@ -53,6 +56,7 @@ plot(awgSignal,'DisplayName','skew 0'); plot(filtered,'DisplayName','filtered'); plot(abs(laserfield),'DisplayName','laser'); plot(abs(att_out),'DisplayName','att_out'); +plot(abs(fib_out),'DisplayName','att_out'); hold off @@ -139,8 +143,4 @@ function yout = applyPulseShaping(xin,fsym,fdac) warning('Check signal length after pulse shaping'); end -end - -function yout = applyBandwidthLimitation(xin) - data_out=ifft(repmat(state.H,1,size(data_in,1)).*fft(data_in.')).'; -end +end \ No newline at end of file