Start Commit

start implementation of class based simulation of a IM/DD communication system. Mostly based on Move-It but cleaned up and with focus on direct detection, however I try to keep the versatility of move-it alive.
This commit is contained in:
Silas Oettinghaus
2023-05-12 15:28:23 +02:00
parent 69df41340f
commit 6d53823466
109 changed files with 1824 additions and 0 deletions

258
Classes/AWG.m Normal file
View File

@@ -0,0 +1,258 @@
classdef AWG
%AWG Summary of this class goes here
% Detailed explanation goes here
properties(Access=public)
kover %oversampling factor e.g. 16
repetitions %repeat the signal to generate a longer sequence?
fdac %needed
normalize %want to normalize at first? either 0 or 1
bit_resolution %bit res. of quantizer (e.g. 5 bit)
dac_min
dac_max
skew_active = 0;
awg_skew ; %skew vector for each output channel
lpf_active = 0;
lpf_type ;
f_cutoff ;
lowpass ;
signal_length;
H_lpf;
end
methods (Access=public)
function obj = AWG(options)
%AWG Construct an instance of this class
% Detailed explanation goes here
arguments
options.preset = [];
options.kover = 16;
options.repetitions = 1;
options.normalize = 1;
options.fdac = 92e9;
options.bit_resolution = 5.5
options.dac_min = -0.5;
options.dac_max = 0.5;
options.skew_active = 0;
options.awg_skew = 0;
options.lpf_active = 0;
options.lpf_type = 0;
options.f_cutoff = 32e9;
options.lowpass = 1;
end
if isempty(options.preset)
obj.skew_active = 0;
obj.lpf_active = 0;
elseif options.preset == "M8196A"
% M8196A (92GBd) https://www.keysight.com/us/en/product/M8196A/92-gsa-s-arbitrary-waveform-generators.html
obj.dac_max = 0.5;
obj.dac_min = -.5;
obj.lowpass = 1; %LP
obj.f_cutoff = 32e9;
elseif options.preset == "M8199B"
%https://www.keysight.com/us/en/assets/3120-1465/data-sheets/M8199A-128-256-GSa-s-Arbitrary-Waveform-Generator.pdf
end
obj.kover = options.kover; %oversampling factor e.g. 16
obj.repetitions = options.repetitions; %repeat the signal to generate a longer sequence?
obj.fdac = options.fdac;
obj.normalize = options.normalize;%want to normalize at first? either 0 or 1
obj.bit_resolution = options.bit_resolution;%bit res. of quantizer (e.g. 5 bit)
obj.dac_min = options.dac_min;
obj.dac_max = options.dac_max;
obj.lpf_active = options.lpf_active;
obj.skew_active = options.skew_active;
obj.awg_skew = options.awg_skew;
end
function data_out = process_channel(obj,data_in)
%METHOD1 Summary of this method goes here
% Detailed explanation goes here
obj.signal_length = length(data_in);
if obj.normalize
% 5.1. Normalize the signal to 1 Vpp and set the amplitude of the signal
data_in = data_in/(max(data_in)-min(data_in));
else
% 5.1. Cut the Signal at -1 and 1 and scale to amplitude
data_in(data_in > 1) = 1;
data_in(data_in < -1) = -1;
end
% 1. Quantize the signal
if obj.bit_resolution>0
data_out = obj.quantization(data_in) ;
else
data_out = data_in;
end
% 2. Sample and hold + repeat (data_out: 1xsignal length)
data_out = repmat(data_out,obj.repetitions,obj.kover);
data_out = reshape(data_out',[],1);
% 3. Add skew
if obj.skew_active
data_out = obj.skew(data_out);
end
% 4. Apply LPF on the signal
if obj.lpf_active
obj.H_lpf = obj.buildFilter(1);
data_out = obj.lpf(data_out) ;
end
end
end
methods (Access=private)
function quant_out = quantization(obj,x_in)
steps = round(2^obj.bit_resolution) + rem(round(2^obj.bit_resolution),2) ;
if isreal(x_in)
% shift signal and clip to quantizer intervall
x_in = min(max(x_in-obj.dac_min,0),obj.dac_max-obj.dac_min);
% quantize signal
x_in = round((steps-1)/(obj.dac_max-obj.dac_min)*x_in);
% scale and shift back
quant_out = x_in*(obj.dac_max-obj.dac_min)/(steps-1) +obj.dac_min;
else
% shift signal and clip to quantizer intervall
x_in = min(max(real(x_in)-obj.dac_min,0),obj.dac_max-obj.dac_min) + ...
1i*min(max(imag(x_in)-obj.dac_min,0),obj.dac_max-obj.dac_min);
% quantize signal and shift back signal
x_in = round((steps-1)/(obj.max-obj.dac_min)*x_in);
% scale and shift back
quant_out = x_in*(obj.dac_max-obj.dac_min)/(steps-1) +(1+1i)*obj.dac_min;
end
end
function skew_out = skew(obj,x_in)
% Exact delay calculation
% Calculate transfer function
fsimu = obj.kover*obj.fdac;
faxis = linspace(fsimu/2, fsimu/2, length(x_in)+1);
faxis = fftshift(faxis(1:end-1));
skew_out = ifft(fft(x_in) .* exp(-1i*2*pi*obj.awg_skew*faxis)' );
skew_out = real(skew_out) ; % get rid of negligible imaginary part
end
function lpf_out = lpf(obj,x_in)
lpf_out = ifft(obj.H_lpf.*fft(x_in));
end
function H = buildFilter(obj,filterType)
filtdegree = 3;
fsimu = obj.kover*obj.fdac;
rp = 0.5; %passband ripple
rs = 0.5; %stopband ripple
switch filterType
case 1
% Bessel filter, impulse invariant transformed
[B, A] = besself(filtdegree, 2*pi*obj.f_cutoff);
[B ,A] = impinvar(B,A,fsimu);
case 2
% Bessel filter, impulse bilinear transformed
[Z, P, K] = besself(filtdegree, 2*pi*obj.f_cutoff);
[Z ,P, K] = bilinear(Z,P,K,fsimu);
[B ,A] = zp2tf(Z ,P ,K);
case 3
% Butterworth filter
if obj.lowpass == 1 %lowpass
[B, A] = butter(filtdegree, obj.f_cutoff/(fsimu/2),'low');
else % highpass
[B, A] = butter(filtdegree, obj.f_cutoff/(fsimu/2),'high');
end
case 4
% Chebyshev 1 filter
[B, A] = cheby1(filtdegree,rp, obj.f_cutoff/(fsimu/2));
case 5
% Chebyshev 2 filter
[B, A] = cheby2(filtdegree,rs, obj.f_cutoff/(fsimu/2));
case 6
% Elliptic filter
[B, A] = ellip(filtdegree,rp,rs,obj.f_cutoff/(fsimu/2));
case 7
% Hamming filter
g=(filtdegree-1)/2;
wc=obj.f_cutoff/(fsimu/2);
B = wc*sinc(wc*(-g:g)).*hamming(filtdegree)';
A=1;
case 8
% Raised Cosine filter
B = firrcos(filtdegree,obj.f_cutoff,para.df,fsimu);
A=1;
case 9
% Sinc filter
g=(filtdegree-1)/2;
wc=obj.f_cutoff/(fsimu/2);
B = wc*sinc(wc*(-g:g));
A=1;
end
H = freqz(B, A, obj.repetitions*obj.kover*obj.signal_length,'whole');
end
end
end

140
Classes/Amplifier.m Normal file
View File

@@ -0,0 +1,140 @@
classdef Amplifier
%UNTITLED Summary of this class goes here
% Detailed explanation goes here
properties
type
amp_mode
amplification_db
nase_mode
noifig
saturation_mode
saturation_power
fsimu
end
methods
function obj = Amplifier(options)
%UNTITLED Construct an instance of this class
% Detailed explanation goes here
arguments
options.type
options.amp_mode
options.amplification_db
options.nase_mode = 0;
options.noifig = 0;
options.saturation_mode = 0;
options.saturation_power = 0;
options.fsimu = []
end
obj.type = options.type;
obj.amp_mode = options.amp_mode;
obj.amplification_db = options.amplification_db;
obj.nase_mode = options.nase_mode;
obj.noifig = options.noifig;
obj.saturation_mode = options.saturation_mode;
obj.saturation_power = options.saturation_power;
obj.fsimu = options.fsimu;
obj.saturation_power=1/1000*10^(obj.saturation_power/10);
end
function [y_out,nase] = process(obj,x_in,optional)
arguments
obj
x_in = [];
optional.nase = 0;
end
%calc gain for output power mode, amp mode and saturated mode
a_lin = obj.calculateGain(x_in);
y_out = a_lin * x_in;
if obj.type == "ideal"
%don't add/ remove noise, but scale it accordingly
if optional.nase ~= 0
nase=state.a^2*optional.nase;
end
elseif obj.type == "edfa"
%calculate ASE-noise
if optional.nase ~= 0
Nase_old = state.gain^2*optional.nase ;
h = Constant.Planck;
c = Constant.LightSpeed;
Nase_new = 0.5* 10^(obj.noifig/10) * h*c/(lambda_T*1e-9) * (a_lin^2-1) ;
nase = Nase_old + Nase_new ;
end
end
pow_in_lin = mean(abs(x_in.^2)) ;
pow_in_dbm = 10*log10(pow_in_lin)+30;
pow_out_lin = mean(abs(y_out.^2)) ;
pow_out_dbm = 10*log10(pow_out_lin)+30;
if obj.amp_mode == "output_power"
seemsright = pow_out_dbm == obj.amplification_db;
else obj.amp_mode == "gain"
seemsright = pow_out_dbm == pow_in_dbm+ obj.amplification_db;
end
if ~seemsright
warning("Amplifier output not correct, please check the reason");
end
end
function a_lin = calculateGain(obj,xin)
if obj.amp_mode == "output_power"
%get linear gain for output power mode
pow_in = mean(abs(xin.^2)) ; % lin input power
pow_out = 10^(obj.amplification_db/10 - 3) ; % dBm to lin
a_lin = sqrt(pow_out/pow_in) ;
elseif obj.amp_mode == "gain"
%get linear gain for classic gain mode
a_lin=10^(obj.amplification_db/20);
end
if obj.saturation_mode
Pin = sum(mean((abs(xin)).^2, 2)) ;
gain_power=fzero(inline(['log(G/' num2str(a_lin^2,'%1.16e') ')/(1-G)-log(2)*' num2str(Pin,'%1.16e') '/' ...
num2str(obj.saturation_power,'%1.16e')],'G'),[1+eps,a_lin^2]);
end
end
end
end

135
Classes/EML.m Normal file
View File

@@ -0,0 +1,135 @@
classdef EML
%EML External Modulated Laser
% Detailed explanation goes here
properties
%class input objms
mode
fsimu %in Hz
lambda %in nm
power %in dBm
linewidth %in Hz
ampl_imbal
pha_imbal
bias
u_pi
%on instance creation
field
noisefactor
phase %eigentlich nur wichtig wenn man mehrere blöcke durchsimuliert...
real_factor
imag_factor
%during simulation
signal_len
end
methods
function obj = EML(options)
%EML Construct an instance of this class
% Detailed explanation goes here
arguments
options.mode;
options.fsimu;
options.lambda;
options.power;
options.linewidth = 0;
options.ampl_imbal = 0;
options.pha_imbal = 0;
options.bias;
options.u_pi;
end
obj.mode = options.mode;
obj.fsimu = options.fsimu;
obj.lambda = options.lambda;
obj.power = options.power;
obj.linewidth = options.linewidth;
obj.ampl_imbal = options.ampl_imbal;
obj.pha_imbal = options.pha_imbal;
obj.bias = options.bias;
obj.u_pi = options.u_pi;
obj.field=sqrt(10^(obj.power/10-3)); %dbm to sqrt(mw)
obj.noisefactor=sqrt(2*pi*obj.linewidth/obj.fsimu);
obj.phase=0;
if obj.mode==4 || obj.mode==5
obj.real_factor=(1+obj.ampl_imbal/2)*exp(1i*obj.pha_imbal/2);
obj.imag_factor=(1-obj.ampl_imbal/2)*exp(1i*(pi/2-obj.pha_imbal/2));
end
end
function [opt_out,obj] = process(obj,elec_in)
%METHOD1 Summary of this method goes here
% Detailed explanation goes here
%needed later
obj.signal_len = length(elec_in);
% generate phase noise of laser
% scale with signal energy in sqrt(mw)
if obj.linewidth ~= 0
ph_noi = obj.createPhaseNoise;
laserfield = obj.field.*exp(1i*ph_noi);
else
laserfield = obj.field;
end
%modulate the laserfield with the electrical signal
opt_out = obj.externalmodulation(laserfield,elec_in);
%remember phase (! you need to receive the altered eml object in you sim program !)
obj.phase = ph_noi(end);
end
function noi = createPhaseNoise(obj)
%create random vector
noi = randn(obj.signal_len,1);
%scale with noisefactor
noi = noi * obj.noisefactor;
%cumsum to accumulate noise over time vector
noi = cumsum(noi);
%add phase from previous block/ loop of simulation
noi = noi + obj.phase;
end
function modulated_laserfield = externalmodulation(obj,laserfield,electrical_in)
switch obj.mode
case 1 % linear IM
modulated_laserfield = laserfield.*sqrt(real(electrical_in + obj.bias)/obj.u_pi);
case 2 % IM with hMZM function
modulated_laserfield = laserfield.*cos(pi/2*(real(electrical_in)+obj.bias)/obj.u_pi);
case 3 % PMsi
modulated_laserfield = laserfield.*exp(1i*pi*real(electrical_in)/obj.u_pi);
case 4 % linear IQ-Modulator
if obj.ampl_imbal ==0 && obj.pha_imbal == 0
% No IQ Imbalance
modulated_laserfield=laserfield.*electrical_in/obj.u_pi;
else
% with IQ Imbalance
modulated_laserfield = laserfield.*(real(electrical_in)*obj.real_factor/obj.u_pi...
+imag(electrical_in)*obj.imag_factor/obj.u_pi);
end
case 5 % IQ-Modulator with cosine function
modulated_laserfield = laserfield/sqrt(2).*(cos(pi/2*(real(electrical_in)+obj.bias)/obj.u_pi)*obj.real_factor...
+cos(pi/2*(imag(electrical_in)+obj.bias)/obj.u_pi)*obj.imag_factor);
end
end
end
end

View File

@@ -0,0 +1,25 @@
classdef Electricalsignal < Signal
%ELECTRICALSIGNAL Summary of this class goes here
% Detailed explanation goes here
properties
end
methods
function obj = Electricalsignal(signal, fsym, fsimu)
%ELECTRICALSIGNAL Construct an instance of this class
% Detailed explanation goes here
obj = obj@Signal(signal, fsym, fsimu);
end
function outputArg = method1(obj,inputArg)
%METHOD1 Summary of this method goes here
% Detailed explanation goes here
outputArg = obj.Property1 + inputArg;
end
end
end

51
Classes/Fiber.m Normal file
View File

@@ -0,0 +1,51 @@
classdef Fiber
%FIBER Summary of this class goes here
% Detailed explanation goes here
properties
fiber_length
alpha
D
Dslope
lambda0
gamma
dphimax
end
methods
function obj = Fiber(options)
%FIBER Construct an instance of this class
% Detailed explanation goes here
arguments
options.fiber_length = 0
options.alpha = 0.2
options.D = 17
options.Dslope = 0.06
options.lambda0 = 1550
options.gamma = 0.0013
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;
end
function opt_out = process(obj,opt_in)
%METHOD1 Summary of this method goes here
% Detailed explanation goes here
%attenuate nase
end
end
end

139
Classes/Filter.m Normal file
View File

@@ -0,0 +1,139 @@
classdef Filter
%FILTER Summary of this class goes here
% Detailed explanation goes here
properties
H
filterType
f_cutoff
signal_length
filtdegree
passband_ripple
stopband_ripple
fdac
end
methods
function obj = Filter(options)
%FILTER Construct an instance of this class
% Detailed explanation goes here
arguments
options.filterType = 1;
options.f_cutoff = 0;
options.fdac = 0;
options.filtdegree = 3;
options.passband_ripple = 0.5;
options.stopband_ripple = 0.5;
end
obj.filterType = options.filterType;
obj.f_cutoff = options.f_cutoff;
obj.filtdegree = options.filtdegree;
obj.passband_ripple = options.passband_ripple;
obj.stopband_ripple = options.stopband_ripple;
obj.fdac = options.fdac;
end
function yout = process(obj,xin)
obj.signal_length = length(xin);
obj.H = obj.buildFilter(obj.filterType);
yout = obj.applyFilter(xin);
end
function y_filtered = applyFilter(obj,xin)
y_filtered = ifft(obj.H.*fft(xin));
end
function H = buildFilter(obj,filterType)
rp = obj.passband_ripple; %passband ripple
rs = obj.stopband_ripple; %stopband ripple
switch filterType
case 1
% Bessel filter, impulse invariant transformed
[B, A] = besself(obj.filtdegree, 2*pi*obj.f_cutoff);
[B ,A] = impinvar(B,A,obj.fdac);
case 2
% Bessel filter, impulse bilinear transformed
[Z, P, K] = besself(obj.filtdegree, 2*pi*obj.f_cutoff);
[Z ,P, K] = bilinear(Z,P,K,obj.fdac);
[B ,A] = zp2tf(Z ,P ,K);
case 3
% Butterworth filter
if obj.lowpass == 1 %lowpass
[B, A] = butter(obj.filtdegree, obj.f_cutoff/(obj.fdac/2),'low');
else % highpass
[B, A] = butter(obj.filtdegree, obj.f_cutoff/(obj.fdac/2),'high');
end
case 4
% Chebyshev 1 filter
[B, A] = cheby1(obj.filtdegree,rp, obj.f_cutoff/(obj.fdac/2));
case 5
% Chebyshev 2 filter
[B, A] = cheby2(obj.filtdegree,rs, obj.f_cutoff/(obj.fdac/2));
case 6
% Elliptic filter
[B, A] = ellip(obj.filtdegree,rp,rs,obj.f_cutoff/(obj.fdac/2));
case 7
% Hamming filter
g=(obj.filtdegree-1)/2;
wc=obj.f_cutoff/(obj.fdac/2);
B = wc*sinc(wc*(-g:g)).*hamming(obj.filtdegree)';
A=1;
case 8
% Raised Cosine filter
B = firrcos(obj.filtdegree,obj.f_cutoff,para.df,obj.fdac);
A=1;
case 9
% Sinc filter
g=(obj.filtdegree-1)/2;
wc=obj.f_cutoff/(obj.fdac/2);
B = wc*sinc(wc*(-g:g));
A=1;
end
H = freqz(B, A, obj.signal_length,'whole');
end
end
end

28
Classes/Opticalsignal.m Normal file
View File

@@ -0,0 +1,28 @@
classdef Opticalsignal < Signal
%OPTICALSIGNAL Summary of this class goes here
% Detailed explanation goes here
properties
nase
lambda_nm
end
methods
function obj = Opticalsignal(signal, fsym, fsimu)
%OPTICALSIGNAL Construct an instance of this class
% Detailed explanation goes here
obj = obj@Signal(signal, fsym, fsimu);
obj.nase = 0;
obj.lambda_nm = [];
end
function outputArg = method1(obj,inputArg)
%METHOD1 Summary of this method goes here
% Detailed explanation goes here
outputArg = obj.Property1 + inputArg;
end
end
end

48
Classes/Signal.m Normal file
View File

@@ -0,0 +1,48 @@
classdef Signal
%SIGNAL Summary of this class goes here
% Detailed explanation goes here
properties
signal
fsym
fsimu
end
methods
function obj = Signal(signal, fsym, fsimu)
%SIGNAL Construct an instance of this class
% Detailed explanation goes here
obj.signal = signal;
obj.fsym = fsym;
obj.fsimu = fsimu;
end
function [e_sig, nase, lambda_nm] = Electricalsignal(obj)
if isa(obj,'Opticalsignal')
%convert to electrical
disp("Convert signal: opt. -> elec. Fetch all properties (e.g. nase)!");
nase = obj.nase;
lambda_nm = obj.lambda_nm;
e_sig = Electricalsignal(obj.signal,obj.fsym, obj.fsimu);
end
end
function o_sig = Opticalsignal(obj)
if isa(obj,'Electricalsignal')
%convert to optical
disp("Convert signal: elec. -> opt.!");
o_sig = Opticalsignal(obj.signal,obj.fsym, obj.fsimu);
end
end
function outputArg = method1(obj,inputArg)
end
end
end

11
Datatypes/emlmodes.m Normal file
View File

@@ -0,0 +1,11 @@
classdef emlmodes < int32
enumeration
im_linear (1)
im_cosinus (2)
pm (3)
iq_linear (4)
iq_cosinus (5)
end
end

15
Datatypes/filtertypes.m Normal file
View File

@@ -0,0 +1,15 @@
classdef filtertypes < int32
enumeration
bessel_inp (1)
bessel_bilin(2)
butterworth (3)
chebyshev1 (4)
chebyshev2 (5)
elliptic (6)
hamming (7)
racos (8)
sinc (9)
end
end

52
Fiber.asv Normal file
View File

@@ -0,0 +1,52 @@
classdef Fiber
%FIBER Summary of this class goes here
% Detailed explanation goes here
properties
fiber_length
alpha
D
Dslope
lambda0
gamma
dphimax
end
methods
function obj = Fiber(options)
%FIBER Construct an instance of this class
% Detailed explanation goes here
arguments
options.fiber_length = 0
options.alpha = 0.2
options.D = 17
options.Dslope = 0.06
options.lambda0 = 1550
options.gamma = 0.0013
options.dphimax = 5e-3
end
obj.fiber_length
obj.alpha
obj.D
Dslope
lambda0
gamma
dphimax
end
function opt_out = process(obj,opt_in)
%METHOD1 Summary of this method goes here
% Detailed explanation goes here
%attenuate nase
end
end
end

31
Signal.asv Normal file
View File

@@ -0,0 +1,31 @@
classdef Signal
%SIGNAL Summary of this class goes here
% Detailed explanation goes here
properties
signal
nase
fsym
fsimu
lambda
end
methods
function obj = Signal(signal, fsym, fsimu)
%SIGNAL Construct an instance of this class
% Detailed explanation goes here
obj.signal = signal;
obj.nase = 0;
obj.fsym = fsym;
obj.fsimu = fsimu;
obj.lambda =
end
function outputArg = method1(obj,inputArg)
%METHOD1 Summary of this method goes here
% Detailed explanation goes here
outputArg = obj.Property1 + inputArg;
end
end
end

146
comm_tb.asv Normal file
View File

@@ -0,0 +1,146 @@
O = 10; %order of prbs
N = 64; %length of prbs
[~,seed] = prbs(O,1); %initialize first seed of prbs
% Modulation
M = 4;
bitpattern = zeros(N,log2(M));
% Symbol Rate
fsym = 56e9;
% DAC Rate
fdac = 120e9;
% Simulation oversampling rate "k";
kover = 16;
% Simulation frequency in "analog domain"
fsimu = kover * fdac ;
%SIMULATE
for i = 1:log2(M)
[bitpattern(:,i),seed] = prbs(O,N,seed);
end
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);
fil = Filter('filtdegree',4,"f_cutoff",50e9,"fdac",fdac,"filterType",filtertypes.bessel_bilin);
filtered = fil.process(awgSignal);
u_pi = 3.5;
vbias = (0.5*u_pi)-u_pi;
eml = EML("mode",emlmodes.im_cosinus,"power",10,"fsimu",fsimu,"lambda",1550,"bias",vbias,"u_pi",u_pi,"linewidth",1000000);
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);
figure;
plot(shapedData);
hold on
plot(awgSignal,'DisplayName','skew 0');
plot(filtered,'DisplayName','filtered');
plot(abs(laserfield),'DisplayName','laser');
plot(abs(att_out),'DisplayName','att_out');
hold off
function pam_sig = pam_mapping(bitpattern, M, unipolar)
switch log2(M)
case 1
% 2-ASK: BPSK / OOK
pam_sig=bitpattern(:,1);
if unipolar==0
pam_sig=2*pam_sig-1;
end
case 2
% 4-ASK:
pam_sig=2*bitpattern(:,1)+(bitpattern(:,1)==bitpattern(:,2));
if unipolar==0
pam_sig=2*pam_sig-3;
end
case 3
% 8-ASK:
x1 = bitpattern(:,1);
x2 = (bitpattern(:,1)==bitpattern(:,3));
x3 = x2~=bitpattern(:,2);
pam_sig = 4*x1 + 2*x2 + x3;
if unipolar==0
pam_sig=2*pam_sig-7;
end
case 4
% 16-ASK:
x1 = bitpattern(:,1);
x2 = (bitpattern(:,1)==bitpattern(:,4));
x3 = x2~=bitpattern(:,3);
x4 = x3~=bitpattern(:,2);
pam_sig = 8*x1 + 4*x2 + 2*x3 + x4;
if unipolar==0
pam_sig=2*pam_sig-15;
end
end
end
function yout = applyPulseShaping(xin,fsym,fdac)
if ~rem(fdac,fsym)
%ist ein Vielfaches
sps = fdac / fsym;
up = sps;
dn = 1;
else
%ist kein Vielfaches
up = fdac / gcd(fdac, fsym);
dn = fsym / gcd(fdac, fsym);
sps= up;
end
%Bau das Filter (hier rrc)
racos_len = 2048;
alpha = 0.1;
h = rcosdesign(alpha,racos_len,sps);
%Apply Filter using Matlab build in fctn.
yout = upfirdn(xin,h,up,dn);
%cut signal, which is longer due to fir filter
st = up/dn*racos_len/2; %we need to cut y_out
en = st + (length(xin)*up/dn) -1;
yout = yout(st:en);
%Check output integrity
if (up/dn * length(xin)) ~= length(yout)
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

146
comm_tb.m Normal file
View File

@@ -0,0 +1,146 @@
O = 10; %order of prbs
N = 64; %length of prbs
[~,seed] = prbs(O,1); %initialize first seed of prbs
% Modulation
M = 4;
bitpattern = zeros(N,log2(M));
% Symbol Rate
fsym = 56e9;
% DAC Rate
fdac = 120e9;
% Simulation oversampling rate "k";
kover = 16;
% Simulation frequency in "analog domain"
fsimu = kover * fdac ;
%SIMULATE
for i = 1:log2(M)
[bitpattern(:,i),seed] = prbs(O,N,seed);
end
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);
fil = Filter('filtdegree',4,"f_cutoff",50e9,"fdac",fdac,"filterType",filtertypes.bessel_bilin);
filtered = fil.process(awgSignal);
u_pi = 3.5;
vbias = (0.5*u_pi)-u_pi;
eml = EML("mode",emlmodes.im_cosinus,"power",10,"fsimu",fsimu,"lambda",1550,"bias",vbias,"u_pi",u_pi,"linewidth",1000000);
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);
figure;
plot(shapedData);
hold on
plot(awgSignal,'DisplayName','skew 0');
plot(filtered,'DisplayName','filtered');
plot(abs(laserfield),'DisplayName','laser');
plot(abs(att_out),'DisplayName','att_out');
hold off
function pam_sig = pam_mapping(bitpattern, M, unipolar)
switch log2(M)
case 1
% 2-ASK: BPSK / OOK
pam_sig=bitpattern(:,1);
if unipolar==0
pam_sig=2*pam_sig-1;
end
case 2
% 4-ASK:
pam_sig=2*bitpattern(:,1)+(bitpattern(:,1)==bitpattern(:,2));
if unipolar==0
pam_sig=2*pam_sig-3;
end
case 3
% 8-ASK:
x1 = bitpattern(:,1);
x2 = (bitpattern(:,1)==bitpattern(:,3));
x3 = x2~=bitpattern(:,2);
pam_sig = 4*x1 + 2*x2 + x3;
if unipolar==0
pam_sig=2*pam_sig-7;
end
case 4
% 16-ASK:
x1 = bitpattern(:,1);
x2 = (bitpattern(:,1)==bitpattern(:,4));
x3 = x2~=bitpattern(:,3);
x4 = x3~=bitpattern(:,2);
pam_sig = 8*x1 + 4*x2 + 2*x3 + x4;
if unipolar==0
pam_sig=2*pam_sig-15;
end
end
end
function yout = applyPulseShaping(xin,fsym,fdac)
if ~rem(fdac,fsym)
%ist ein Vielfaches
sps = fdac / fsym;
up = sps;
dn = 1;
else
%ist kein Vielfaches
up = fdac / gcd(fdac, fsym);
dn = fsym / gcd(fdac, fsym);
sps= up;
end
%Bau das Filter (hier rrc)
racos_len = 2048;
alpha = 0.1;
h = rcosdesign(alpha,racos_len,sps);
%Apply Filter using Matlab build in fctn.
yout = upfirdn(xin,h,up,dn);
%cut signal, which is longer due to fir filter
st = up/dn*racos_len/2; %we need to cut y_out
en = st + (length(xin)*up/dn) -1;
yout = yout(st:en);
%Check output integrity
if (up/dn * length(xin)) ~= length(yout)
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

View File

@@ -0,0 +1,6 @@
function [v,u]=Boltzmann()
v=1;
u=9.1e-7;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=ConductanceQuantum()
v=1/pi;
u=1e-10;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=Coulomb()
v=1;
u=0;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=ElectronMass()
v=1;
u=4.4e-8;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=ElementaryCharge()
v=1;
u=2.2e-8;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=FineStructure()
v=7.2973525698e-3;
u=2.4e-12;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=FluxQuantum()
v=pi;
u=6.9e-8;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=LightSpeed()
v=137.03599907;
u=4.4e-8;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=Planck()
v=2*pi;
u=7.5e-8;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=ProtonMass()
v=1836.15267;
u=1.6e-4;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=ReducedPlanck()
v=1;
u=1.2e-8;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=Boltzmann()
v=1;
u=9.1e-7;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=Coulomb()
v=1;
u=0;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=Gravitational()
v=1;
u=4.7e-5;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=LightSpeed()
v=1;
u=0;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=Planck()
v=2*pi;
u=7.5e-8;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=ReducedPlanck()
v=1;
u=1.2e-8;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=Boltzmann()
v=1;
u=9.1e-7;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=LightSpeed()
v=1;
u=0;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=Planck()
v=2*pi;
u=7.5e-8;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=ProtonMass()
v=1;
u=4.4e-8;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=ReducedPlanck()
v=1;
u=1.2e-8;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=Boltzmann()
v=1;
u=9.1e-7;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=ConductanceQuantum()
v=2/pi;
u=2e-10;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=Coulomb()
v=1;
u=0;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=ElectronMass()
v=0.5;
u=2.2e-8;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=ElementaryCharge()
v=sqrt(2);
u=3.1e-8;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=FineStructure()
v=7.2973525698e-3;
u=2.4e-12;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=FluxQuantum()
v=pi/sqrt(2);
u=4.9e-8;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=LightSpeed()
v=274.07199814;
u=8.8e-8;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=Planck()
v=2*pi;
u=7.5e-8;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=ProtonMass()
v=918.076335;
u=8.1e-5;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=ReducedPlanck()
v=1;
u=1.2e-8;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=Boltzmann()
v=1;
u=9.1e-7;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=Coulomb()
v=1;
u=0;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=ElementaryCharge()
v=1;
u=2.2e-8;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=Gravitational()
v=1;
u=4.7e-5;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=LightSpeed()
v=1;
u=0;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=AtomicMass()
v=1.660538921e-27;
u=7.3e-35;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=AvogadroNumber()
v=6.02214129e23;
u=2.7e16;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=BohrMagneton()
v=9.27400968e-24;
u=2e-31;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=BohrRadius()
v=5.2917721092e-11;
u=1.7e-20;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=Boltzmann()
v=1.3806488e-23;
u=1.3e-29;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=ConductanceQuantum()
v=7.7480917346e-5;
u=2.5e-14;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=Coulomb()
v=8.987551787;
u=0;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=ElectronMass()
v=9.10938291e-31;
u=4e-38;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=ElementaryCharge()
v=1.602176565e-19;
u=3.5e-27;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=Faraday()
v=96485.3365;
u=2.1e-3;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=FermiCoupling()
v=1.166364e-5;
u=5e-11;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=FineStructure()
v=7.2973525698e-3;
u=2.4e-12;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=FluxQuantum()
v=2.067833758e-15;
u=4.6e-23;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=Gravitational()
v=6.67408e-11;
u=3.1e-15;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=HartreeEnergy()
v=4.35974434e-18;
u=1.9e-25;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=LightSpeed()
v=299792458;
u=0;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=Loschmidt()
v=2.6867805e25;
u=2.4e19;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=NuclearMagneton()
v=5.05078353e-27;
u=1.1e-34;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=Planck()
v=6.626070040e-34;
u=8.1e-42;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=ProtonMass()
v=1.672621777e-27;
u=7.4e-35;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=ReducedPlanck()
v=1.054571800e-34;
u=1.3e-42;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=Rydberg()
v=10973731.568539;
u=5.5e-5;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=StefanBoltzmann()
v=5.670373e-8;
u=2.1e-13;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=VacuumImpedance()
v=376.730313461;
u=0;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=VacuumPermeability()
v=1.256637061e-6;
u=0;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=VacuumPermittivity()
v=8.854187817e-12;
u=0;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=VonKlitzing()
v=25812.8074434;
u=8.4e-6;
end

View File

@@ -0,0 +1,6 @@
function [v,u]=WienDisplacement()
v=2.8977721e-3;
u=2.6e-9;
end

View File

@@ -0,0 +1,5 @@
function [v,u]=AtomicMass()
[v,u]=Constant.(UnitSystem()).AtomicMass;
end

View File

@@ -0,0 +1,5 @@
function [v,u]=AvogadroNumber()
[v,u]=Constant.(UnitSystem()).AvogadroNumber;
end

View File

@@ -0,0 +1,5 @@
function [v,u]=BohrMagneton()
[v,u]=Constant.(UnitSystem()).BohrMagneton;
end

View File

@@ -0,0 +1,5 @@
function [v,u]=BohrRadius()
[v,u]=Constant.(UnitSystem()).BohrRadius;
end

View File

@@ -0,0 +1,5 @@
function [v,u]=Boltzmann()
[v,u]=Constant.(UnitSystem()).Boltzmann;
end

View File

@@ -0,0 +1,5 @@
function [v,u]=ConductanceQuantum()
[v,u]=Constant.(UnitSystem()).ConductanceQuantum;
end

View File

@@ -0,0 +1,5 @@
function [v,u]=Coulomb()
[v,u]=Constant.(UnitSystem()).Coulomb;
end

View File

@@ -0,0 +1,5 @@
function [v,u]=ElectronMass()
[v,u]=Constant.(UnitSystem()).ElectronMass;
end

View File

@@ -0,0 +1,5 @@
function [v,u]=ElementaryCharge()
[v,u]=Constant.(UnitSystem()).ElementaryCharge;
end

View File

@@ -0,0 +1,5 @@
function [v,u]=Faraday()
[v,u]=Constant.(UnitSystem()).Faraday;
end

View File

@@ -0,0 +1,5 @@
function [v,u]=FermiCoupling()
[v,u]=Constant.(UnitSystem()).FermiCoupling;
end

View File

@@ -0,0 +1,5 @@
function [v,u]=FineStructure()
[v,u]=Constant.(UnitSystem()).FineStructure;
end

View File

@@ -0,0 +1,5 @@
function [v,u]=FluxQuantum()
[v,u]=Constant.(UnitSystem()).FluxQuantum;
end

View File

@@ -0,0 +1,5 @@
function [v,u]=Gravitational()
[v,u]=Constant.(UnitSystem()).Gravitational;
end

View File

@@ -0,0 +1,5 @@
function [v,u]=HartreeEnergy()
[v,u]=Constant.(UnitSystem()).HartreeEnergy;
end

View File

@@ -0,0 +1,5 @@
function [v,u]=LightSpeed()
[v,u]=Constant.(UnitSystem()).LightSpeed;
end

View File

@@ -0,0 +1,5 @@
function [v,u]=Loschmidt()
[v,u]=Constant.(UnitSystem()).LoschmidtConstant;
end

View File

@@ -0,0 +1,5 @@
function [v,u]=NuclearMagneton()
[v,u]=Constant.(UnitSystem()).NuclearMagneton;
end

View File

@@ -0,0 +1,5 @@
function [v,u]=Planck()
[v,u]=Constant.(UnitSystem()).Planck;
end

View File

@@ -0,0 +1,5 @@
function [v,u]=ProtonMass()
[v,u]=Constant.(UnitSystem()).ProtonMass;
end

Some files were not shown because too many files have changed in this diff Show More