WIP - implementation and debugging of 224 Gbd PAM4 system

Classes/AWG.m

@@ -18,7 +18,6 @@ classdef AWG
         lpf_active = 0;
         lpf_type ;
         f_cutoff ;
-        lowpass ;
 
         signal_length;
         H_lpf;
@@ -45,35 +44,30 @@ classdef AWG
                 options.lpf_active      = 0;
                 options.lpf_type        = 0;
                 options.f_cutoff        = 32e9;
-                options.lowpass         = 1;
             end
 
-            if isempty(options.preset)
+            fn = fieldnames(options);
-                obj.skew_active = 0;
+            for n = 1:numel(fn)
-                obj.lpf_active = 0;
+                try
+                    obj.(fn{n}) = options.(fn{n});
+                end
+            end
 
-            elseif options.preset == "M8196A"
+            if 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        = 50e9;
             elseif options.preset == "M8199B"
                 %https://www.keysight.com/us/en/assets/3120-1465/data-sheets/M8199A-128-256-GSa-s-Arbitrary-Waveform-Generator.pdf
+                % obj.fdac            = 256e9;
+                obj.dac_max         = 0.5;
+                obj.dac_min         = -.5;
+                obj.f_cutoff        = 80e9;
+
             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;
 
 
 
@@ -81,22 +75,22 @@ classdef AWG
 
         function signalclass_out = process(obj,signalclass_in)
 
-             % actual processing of the signal (steps 1. - 3.)
+            % actual processing of the signal (steps 1. - 3.)
             signalclass_in.signal = obj.process_(signalclass_in.signal);
 
             % 4. Apply LPF on the signal
             if obj.lpf_active
-                lpf = Filter('filtdegree',4,"f_cutoff",obj.f_cutoff,"fsamp",obj.kover*obj.fdac,"filterType",filtertypes.bessel_inp);
+                lpf = Filter('filtdegree',5,"f_cutoff",obj.f_cutoff,"fsamp",obj.kover*obj.fdac,"filterType",obj.lpf_type);
                 signalclass_in = lpf.process(signalclass_in);
             end
 
-             % cast the inform. signal to electrical signal 
+            % cast the inform. signal to electrical signal
             signalclass_in = Electricalsignal(signalclass_in,"fs",obj.fdac*obj.kover,"logbook",signalclass_in.logbook);
 
-             % append to logbook
+            % append to logbook
             signalclass_in  = signalclass_in.logbookentry();
 
-             % write to output
+            % write to output
             signalclass_out = signalclass_in;
 
         end
@@ -116,15 +110,14 @@ classdef AWG
             obj.signal_length = length(data_in);
 
             if obj.normalize
-                % 5.1. Normalize the signal to 1 Vpp and set the amplitude of the signal
+                % 0a 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
+                % 0b 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
                 elec_out = obj.quantization(data_in) ;
@@ -141,8 +134,6 @@ classdef AWG
                 elec_out = obj.skew(elec_out);
             end
 
-
-
         end
 
     end
@@ -194,92 +185,5 @@ classdef AWG
 
         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

Classes/Amplifier.m

@@ -3,18 +3,14 @@ classdef Amplifier
     %   Detailed explanation goes here
 
     properties
-        type 
         amp_mode
-        amplification_db
+        gain_mode
-
         nase_mode
+
+        amplification_db       
         noifig
 
-        saturation_mode
-        saturation_power
-
         fsimu
-
     end
 
     methods
@@ -22,134 +18,145 @@ classdef Amplifier
             %UNTITLED Construct an instance of this class
             %   Detailed explanation goes here
             arguments
-                options.type
+                options.amp_mode amp_mode   = amp_mode.ideal_no_noise
-                options.amp_mode
+                options.gain_mode gain_mode = gain_mode.output_power
-                options.amplification_db
+                options.nase_mode nase_mode = nase_mode.pass_ase
 
-                options.nase_mode = 0;
+                options.amplification_db double = 0
-                options.noifig = 0;
                 
-                options.saturation_mode = 0;
+                options.noifig double = 0;
-                options.saturation_power = 0;
-
-                options.fsimu = []
             end
 
-            obj.type = options.type;
+            fn = fieldnames(options);
-            obj.amp_mode = options.amp_mode;
+            for n = 1:numel(fn)
-            obj.amplification_db = options.amplification_db;
+                obj.(fn{n}) = options.(fn{n});
-
+            end
-            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 signalclass_out = process(obj,signalclass_in)
+        function signalclass_out = process(obj,signalclass_in)
 
-             % actual processing of the signal (steps 1. - 3.)
+            % actual processing of the signal
-            signalclass_in.signal = obj.process_(signalclass_in.signal); 
+            signalclass_in = obj.process_(signalclass_in);
 
-             % append to logbook
+            % append to logbook
-             lbdesc = ['Amp '];
+            lbdesc = ['Amp '];
             signalclass_in  = signalclass_in.logbookentry(lbdesc);
 
-             % write to output
+            % write to output
             signalclass_out = signalclass_in;
 
         end
 
+        function [X_out] = process_(obj,X_in)
 
-        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);
+            a_lin = obj.calculateGain(X_in.signal);
 
-            y_out = a_lin * x_in;
+            %apply amplification
+            X_in.signal = a_lin * X_in.signal;
             
-            if obj.type == "ideal"
+            % now, optical ase noise:
+            if class(X_in) == "Opticalsignal"
+
+                if obj.amp_mode == amp_mode.ideal_no_noise
+
+                    %don't add/ remove noise, but scale it accordingly
+                    X_in.nase = obj.onlyAmplifyAse(X_in.nase, a_lin);
+
+                elseif obj.amp_mode == amp_mode.edfa_increase_nase
+
+                    %calculate ASE-noise
+                    X_in.nase = obj.increaseAse(X_in.nase, a_lin, obj.noifig , X_in.lambda );
+
+                elseif obj.amp_mode == amp_mode.edfa_replace_nase
+
+                    %calculate ASE-noise
+                    X_in.nase = obj.replaceAse(X_in.nase, a_lin, obj.noifig , X_in.lambda );
                     
-                %don't add/ remove noise, but scale it accordingly
-                if optional.nase ~= 0
-                    nase=state.a^2*optional.nase;
                 end
 
-            elseif obj.type == "edfa"
+                if obj.nase_mode == nase_mode.generate_ase
                     
-                %calculate ASE-noise
+                    nase        = X_in.nase;
-                if optional.nase ~= 0
+                    fs          = X_in.fs;
-                    Nase_old = state.gain^2*optional.nase ;
+                    dimension   = size(X_in.signal);
 
-                    h = Constant.Planck;
+                    nase_numeric = obj.generateAseNoise(nase, fs, dimension);
-                    c = Constant.LightSpeed;
 
-                    Nase_new = 0.5* 10^(obj.noifig/10) * h*c/(lambda_T*1e-9) * (a_lin^2-1) ;
+                    [X_in.signal, osnr] = obj.applyAseToSignal(X_in.signal, nase_numeric);
+
+                    X_in.nase   = 0;
+
+                elseif obj.nase_mode == nase_mode.pass_ase
+                    
+                    X_in.nase = X_in.nase;
 
-                    nase = Nase_old + Nase_new ; 
                 end
 
             end
 
-            pow_in_lin =  mean(abs(x_in.^2)) ; 
+            X_out = X_in;
-            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;
-            elseif 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"
+            if obj.gain_mode == gain_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"
+            elseif obj.gain_mode == gain_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
         
+        function nase_amped = onlyAmplifyAse(~,nase_old, a_lin)
+            nase_amped = a_lin^2 * nase_old;
+        end
 
+        function nase_new = calculateAseFromThisAmp(~, a_lin, noisefig, lambda)
+            h = Constant.Planck;
+            c = Constant.LightSpeed;
+            nase_new = 0.5* 10^(noisefig/10) * h*c/lambda * (a_lin^2-1) ;
+        end
 
+        function ase_increased = increaseAse(obj,nase_old, a_lin, noisefig, lambda)           
+            
+            nase_fromThisAmp = obj.calculateAseFromThisAmp(a_lin, noisefig, lambda);
+            
+            nase_old = a_lin^2 * nase_old;
+
+            ase_increased = nase_old + nase_fromThisAmp;
+        end
+
+        function nase_new = replaceAse(~,a_lin, noisefig, lambda_nm)
+
+            nase_fromThisAmp = obj.calculateAseFromThisAmp(a_lin, noisefig, lambda_nm);
+            
+            nase_new = nase_fromThisAmp;
+        end
+
+        function nase_numeric = generateAseNoise(~, nase, fs, dimension)
+            nase_numeric = (randn(dimension) + 1i*randn(dimension))*sqrt(nase/2*fs) ;   
+        end
+
+        function [noisy_sig, osnr] = applyAseToSignal(~, opticalsignal, nase_numeric)
+            noisy_sig = opticalsignal + nase_numeric ;
+            osnr = pow2db(mean(abs(opticalsignal.^2))/mean(abs(nase_numeric.^2)));
+            disp(osnr);
+        end
     end
 end
 

Classes/Filter.m

@@ -6,13 +6,15 @@ classdef Filter
         H
         filterType
         f_cutoff
-
+        lowpass
 
         signal_length
         filtdegree
         passband_ripple
         stopband_ripple
         fsamp
+
+        w
     end
 
     methods
@@ -20,36 +22,34 @@ classdef Filter
             %FILTER Construct an instance of this class
             %   Detailed explanation goes here
             arguments
-                options.filterType = 1;
+                options.filterType filtertypes = filtertypes.bessel_inp ;
                 options.f_cutoff   = 0;
                 options.fsamp      = 0;
                 options.filtdegree = 3;
                 options.passband_ripple = 0.5;
                 options.stopband_ripple = 0.5;
+                options.lowpass    = 1;
             end
 
 
-            obj.filterType = options.filterType;
+            fn = fieldnames(options);
-
+            for n = 1:numel(fn)
-            obj.f_cutoff        = options.f_cutoff;
+                obj.(fn{n}) = options.(fn{n});
-            obj.filtdegree      = options.filtdegree;
+            end
-            obj.passband_ripple = options.passband_ripple;
-            obj.stopband_ripple = options.stopband_ripple;
-            obj.fsamp           = options.fsamp;
 
 
         end
 
-         function signalclass_out = process(obj,signalclass_in)
+        function signalclass_out = process(obj,signalclass_in)
 
-             % actual processing of the signal 
+            % actual processing of the signal
             signalclass_in.signal = obj.process_(signalclass_in.signal);
 
-             % append to logbook
+            % append to logbook
             filterdesc = [num2str(obj.filtdegree),'. order ',char(obj.filterType),' filter with f_cutoff at ', num2str(obj.f_cutoff*1e-9), ' GHz.'];
             signalclass_in  = signalclass_in.logbookentry(filterdesc);
 
-             % write to output
+            % write to output
             signalclass_out = signalclass_in;
 
         end
@@ -58,7 +58,7 @@ classdef Filter
 
             obj.signal_length = length(xin);
 
-            obj.H = obj.buildFilter(obj.filterType);
+            [obj.H,obj.w] = obj.buildFilter(obj.filterType);
 
             yout = obj.applyFilter(xin);
 
@@ -70,8 +70,8 @@ classdef Filter
 
         end
 
-        function H = buildFilter(obj,filterType)
+        function [H,w] = buildFilter(obj,filterType)
-
+            w = [];
             rp = obj.passband_ripple; %passband ripple
             rs = obj.stopband_ripple; %stopband ripple
 
@@ -143,9 +143,22 @@ classdef Filter
                     B = wc*sinc(wc*(-g:g));
                     A=1;
 
+                case 10
+
+                    % Gaussian Filter
+                    %check if order ist multiple of 1/2
+                    blocklen=obj.signal_length;
+            
+                    faxis=linspace(-obj.fsamp/2,obj.fsamp/2,blocklen+1)';%generates arow vector faxis of blocklen+1 points linearly spaced between and including -para.fs/2 and para.fs/2
+                    faxis=ifftshift(faxis(1:end-1));
+                    
+                    H=exp(-((faxis)/(obj.f_cutoff*2)).^(2*obj.filtdegree)*log(2)*2^(2*obj.filtdegree-1));
             end
 
-            H = freqz(B, A, obj.signal_length,'whole');
+            % Build Filter from coefficients
+            if filterType ~= 10
+                [H,w] = freqz(B, A, obj.signal_length,'whole');
+            end
 
         end
     end

Classes/Scope.m

@@ -71,7 +71,7 @@ classdef Scope
         function signalclass_out = process(obj,signalclass_in)
 
             % apply LPF
-            lpf = Filter('filtdegree',4,"f_cutoff",obj.lpf_bw,"fsamp",obj.fsimu,"filterType",obj.filtertype);
+            lpf = Filter('filtdegree',3,"f_cutoff",obj.lpf_bw,"fsamp",obj.fsimu,"filterType",obj.filtertype);
             signalclass_in = lpf.process(signalclass_in);
 
             % actual processing of the signal

Classes/Signal.m

@@ -91,7 +91,6 @@ classdef Signal
             if isa(obj,'Electricalsignal')
 
                 %convert to optical
-                disp("Convert signal: elec. -> opt.!");
                 
                 o_sig = Opticalsignal(obj.signal,"fs",obj.fs,"lambda",options.lambda,"logbook",obj.logbook,"nase",options.nase);
  

Datatypes/amp_mode.m

@@ -0,0 +1,10 @@
+classdef amp_mode < int32
+
+   enumeration
+      ideal_no_noise     (1)
+      edfa_increase_nase (2)
+      edfa_replace_nase  (3)
+      % edfa_set_osnr      (4) TODO: implement mode to achieve desired OSNR
+   end
+   
+end

Datatypes/filtertypes.m

@@ -10,6 +10,7 @@ classdef filtertypes < int32
       hamming     (7)
       racos       (8)
       sinc        (9)
+      gaussian    (10)
    end
 
 end

Datatypes/gain_mode.m

@@ -0,0 +1,8 @@
+classdef gain_mode < int32
+
+   enumeration
+      gain     (1)
+      output_power (2)
+   end
+   
+end

Datatypes/nase_mode.m

@@ -0,0 +1,8 @@
+classdef nase_mode < int32
+
+   enumeration
+      generate_ase (1)
+      pass_ase (2)
+   end
+   
+end

Functions/functionSignatures.json

@@ -0,0 +1,14 @@
+{
+  "_schemaVersion": "1.0.0",
+  "myFunc":
+  {
+     "inputs":
+     [
+        {"name":"in1", "kind":"required", "type":["numeric"], "purpose":"ID of item"},
+        {"name":"in2", "kind":"required", "type":["numeric"], "purpose":"# Items"},
+        {"name":"in3", "kind":"ordered", "type":["numeric"], "purpose":"Input Value"},
+        {"name":"Name1", "kind":"namevalue", "type":["logical","scalar"],"purpose":"Option"},
+        {"name":"Name2", "kind":"namevalue", "type":["char", "choices={'Default','Choice1','Choice2'}"]}
+     ]
+  }
+}

Functions/myFunc.m

@@ -0,0 +1,34 @@
+% myFunc  Example function
+% This function is called with any of these syntaxes:
+%
+%   myFunc(in1, in2) accepts 2 required arguments. 
+%   myFunc(in1, in2, in3) also accepts an optional 3rd argument. 
+%   myFunc(___, NAME, VALUE) accepts one or more of the following name-value pair 
+%       arguments. This syntax can be used in any of the previous syntaxes.
+%           * 'NAME1' with logical value
+%           * 'NAME2' with 'Default', 'Choice1', or 'Choice2'
+function myFunc(reqA,reqB,varargin)
+    % Initialize default values
+    NV1 = true;
+    NV2 = 'Default';
+    posA = [];
+    
+    if nargin > 3
+        if rem(nargin,2)
+            posA = varargin{1};
+            V = varargin(2:end);
+        else
+            V = varargin;
+        end
+        for n = 1:2:size(V,2)
+            switch V{n}
+                case 'Name1'
+                    NV1 = V{n+1};
+                case 'Name2'
+                    NV2 = V{n+1}
+                otherwise
+                    error('Error.')
+            end
+        end
+    end
+end

main_simulation.m

@@ -1,6 +1,6 @@
-clear all
+%clear all
 
-O = 15; %order of prbs
+O = 17; %order of prbs
 N = 2^(O-1); %length of prbs
 [~,seed] = prbs(O,1); %initialize first seed of prbs
 
@@ -9,7 +9,7 @@ M = 4; %PAM-M
 bitpattern = zeros(N,log2(M));
 
 % Symbol Rate
-fsym = 56e9;
+fsym = 92e9;
 % DAC Rate
 fdac = 120e9;
 % Simulation oversampling rate "k";
@@ -22,25 +22,33 @@ fsimu = kover * fdac ;
 %CONSTRUCTION
 pam_mapper = PAMmapper(M,0);
 
-awg = AWG('preset','M8196A','fdac',fdac,'kover',kover,'lpf_active',0,'f_cutoff',80e9,'lpf_type',filtertypes.bessel_bilin);
+awg = AWG('preset','M8199B','fdac',fdac,'kover',kover,'lpf_active',1,'f_cutoff',80e9,'lpf_type',filtertypes.gaussian);
 
-%fil_tx = Filter('filtdegree',1,"f_cutoff",80e9,"fsamp",fdac,"filterType",filtertypes.bessel_bilin);
+
-fil_tx = Filter('filtdegree',1,"f_cutoff",50e9,"fsamp",fdac,"filterType",filtertypes.bessel_bilin);
+fil_tx = Filter('filtdegree',1,"f_cutoff",50e9,"fsamp",fdac,"filterType",filtertypes.bessel_inp);
 
 u_pi = 3.5;
 vbias = (0.5*u_pi)-u_pi;
-extmodlaser = EML("mode",emlmodes.im_cosinus,"power",10,"fsimu",fsimu,"lambda",1550,"bias",vbias,"u_pi",u_pi,"linewidth",1000);
+extmodlaser = EML("mode",emlmodes.im_cosinus,"power",5,"fsimu",fsimu,"lambda",1550,"bias",vbias,"u_pi",u_pi,"linewidth",0);
 
-amp = Amplifier("amplification_db",10,"amp_mode","gain","type","ideal","saturation_mode",0,'saturation_power',10);
+amp = Amplifier("amp_mode","ideal_no_noise","amplification_db",0,"gain_mode","output_power");
+
+fib = Fiber("fsimu",fdac*kover,"fiber_length",0,"alpha",0.2,"D",17,"lambda0",1550);
+
+%noise loading 
+zeroDB_attenuator = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",0);
+
+edfa_             = Amplifier("amp_mode","edfa_increase_nase","nase_mode","generate_ase","amplification_db",0,"noifig",33,"gain_mode","output_power");
+
+rop_amplifier     = Amplifier("amp_mode","ideal_no_noise","amplification_db",0,"gain_mode","output_power","nase_mode","pass_ase");
 
-fib = Fiber("fsimu",fdac*kover,"fiber_length",2,"alpha",0.2,"D",17,"lambda0",1550);
 
 phdiode = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20);
 
-fil_diode = Filter('filtdegree',4,"f_cutoff",70e9,"fsamp",fdac,"filterType",filtertypes.bessel_bilin);
+fil_diode = Filter('filtdegree',1,"f_cutoff",70e9,"fsamp",fdac,"filterType",filtertypes.bessel_inp);
 
 scp = Scope("fsimu",fdac*kover,"fadc",fadc,...
-    "delay",0,"fixed_delay",0,"lpf_bw",120e9,"filtertype",filtertypes.bessel_inp,...
+    "delay",0,"fixed_delay",0,"lpf_bw",113e9,"filtertype",filtertypes.butterworth,...
     "samplingdelay",0,"rand_samplingdelay",0,"freq_offset",0,"samp_jitter",0,...
     "adcresolution",6,"quantbuffer",0.1);
 
@@ -48,24 +56,24 @@ eq = EQ("K",2,"plottrain",0,"plotfinal",1,"training_length",2048,"Ne",[25,5,5],"
 
 %SIMULATE
 
-
 % INFORMATION SIGNAL
 for i = 1:log2(M)
     [bitpattern(:,i),seed] = prbs(O,N,seed);
 end
 
-X = Informationsignal(bitpattern);
+bits = Informationsignal(bitpattern);
 
-PAMSIG = pam_mapper.map(X);
+mod_out = pam_mapper.map(bits);
-
+reference = mod_out;
-X = PAMSIG;
+mod_out.signal  = applyPulseShaping(mod_out.signal,fsym,fdac);
-
-X.signal  = applyPulseShaping(PAMSIG.signal,fsym,fdac);
 
 % ELECTRICAL DOMAIN
-X  = awg.process(X);
+X  = awg.process(mod_out);
 
-X  = fil_tx.process(X);
+% X  = fil_tx.process(X);
+%X  = fil_tx.process(X);
+
+X = X.normalize;
 
 % OPTICAL DOMAIN
 X  = extmodlaser.process(X);
@@ -74,6 +82,10 @@ X  = amp.process(X);
 
 X  = fib.process(X);
 
+% X = zeroDB_attenuator.process(X);
+% X = edfa_.process(X);
+% X = rop_amplifier.process(X);
+
 X  = phdiode.process(X);
 
 X  = fil_diode.process(X);
@@ -84,15 +96,25 @@ X  = scp.process(X);
 X  = X.resample("fs_out",2*fsym,"fs_in",fadc);
 
 % INFORMATION SIGNAL
-X  = eq.process(X,PAMSIG);
+X  = X.normalize;
+X  = eq.process(X,reference);
+rx_series = X.signal;
 
 X  = pam_mapper.demap(X);
 
+
 % BER
-[bits,errors,BER] = calc_ber(X.signal(:,10000:end-21),bitpattern(10000:end-20,:)',0);
+[bits,errors,BER] = calc_ber(X.signal(:,10000:end-20),bitpattern(10000:end-20,:)',0);
 
 disp(['BER: ', sprintf('%2E',BER)]);
 
+
+figure()
+hold on
+plot(reference.signal);
+plot(rx_series);
+hold off
+
 disp(X.logbook);
 
 function yout = applyPulseShaping(xin,fsym,fdac)