WIP

Classes/AWG.m

@@ -57,7 +57,7 @@ classdef AWG
                 obj.dac_max         = 0.5;
                 obj.dac_min         = -.5;
                 obj.lowpass         = 1; %LP
-                obj.f_cutoff        = 32e9;
+                obj.f_cutoff        = 92e9;
             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
@@ -112,8 +112,10 @@ classdef AWG
 
             % 4. Apply LPF on the signal
             if obj.lpf_active
-                obj.H_lpf = obj.buildFilter(1);
+                lpf = Filter('filtdegree',4,"f_cutoff",obj.f_cutoff,"fsamp",obj.kover*obj.fdac,"filterType",filtertypes.bessel_inp);
-                data_out = obj.lpf(data_out) ;
+                data_out = lpf.process(data_out);
+%                 obj.H_lpf = obj.buildFilter(1);
+%                 data_out = obj.lpf(data_out) ;
             end
 
         end

Classes/Filter.m

@@ -12,7 +12,7 @@ classdef Filter
         filtdegree
         passband_ripple 
         stopband_ripple
-        fdac
+        fsamp
     end
 
     methods
@@ -22,7 +22,7 @@ classdef Filter
             arguments
                 options.filterType = 1;
                 options.f_cutoff   = 0;
-                options.fdac       = 0;
+                options.fsamp      = 0;
                 options.filtdegree = 3;
                 options.passband_ripple = 0.5;
                 options.stopband_ripple = 0.5;
@@ -31,12 +31,12 @@ classdef Filter
 
             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;
+            obj.fsamp           = options.fsamp;
+            
 
         end
 
@@ -67,49 +67,49 @@ classdef Filter
                     % Bessel filter, impulse invariant transformed
 
                     [B, A] = besself(obj.filtdegree, 2*pi*obj.f_cutoff);
-                    [B ,A] = impinvar(B,A,obj.fdac);
+                    [B ,A] = impinvar(B,A,obj.fsamp);
 
                 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);
+                    [Z ,P, K] = bilinear(Z,P,K,obj.fsamp);
                     [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');
+                        [B, A] = butter(obj.filtdegree, obj.f_cutoff/(obj.fsamp/2),'low');
                     else % highpass
-                        [B, A] = butter(obj.filtdegree, obj.f_cutoff/(obj.fdac/2),'high');
+                        [B, A] = butter(obj.filtdegree, obj.f_cutoff/(obj.fsamp/2),'high');
                     end
 
                 case 4
 
                     % Chebyshev 1 filter
 
-                    [B, A] = cheby1(obj.filtdegree,rp, obj.f_cutoff/(obj.fdac/2));
+                    [B, A] = cheby1(obj.filtdegree,rp, obj.f_cutoff/(obj.fsamp/2));
 
                 case 5
 
                     % Chebyshev 2 filter
 
-                    [B, A] = cheby2(obj.filtdegree,rs, obj.f_cutoff/(obj.fdac/2));
+                    [B, A] = cheby2(obj.filtdegree,rs, obj.f_cutoff/(obj.fsamp/2));
 
                 case 6
 
                     % Elliptic filter
 
-                    [B, A] = ellip(obj.filtdegree,rp,rs,obj.f_cutoff/(obj.fdac/2));
+                    [B, A] = ellip(obj.filtdegree,rp,rs,obj.f_cutoff/(obj.fsamp/2));
 
                 case 7
 
                     % Hamming filter
 
                     g=(obj.filtdegree-1)/2;
-                    wc=obj.f_cutoff/(obj.fdac/2);
+                    wc=obj.f_cutoff/(obj.fsamp/2);
                     B = wc*sinc(wc*(-g:g)).*hamming(obj.filtdegree)';
                     A=1;
 
@@ -117,7 +117,7 @@ classdef Filter
 
                     % Raised Cosine filter
 
-                    B = firrcos(obj.filtdegree,obj.f_cutoff,para.df,obj.fdac);
+                    B = firrcos(obj.filtdegree,obj.f_cutoff,para.df,obj.fsamp);
                     A=1;
 
                 case 9
@@ -125,7 +125,7 @@ classdef Filter
                     % Sinc filter
 
                     g=(obj.filtdegree-1)/2;
-                    wc=obj.f_cutoff/(obj.fdac/2);
+                    wc=obj.f_cutoff/(obj.fsamp/2);
                     B = wc*sinc(wc*(-g:g));
                     A=1;
 

Classes/Photodiode.m

@@ -37,10 +37,10 @@ classdef Photodiode
             R = 50; %resistance of phdiode (50ohm is typical value)
 
             % Magnitude squared detection
-            yout = sum( abs(xin) .^2*obj.responsivity,1) ;
+            yout = sum( abs(xin) .^2*obj.responsivity, 2 ) ;
 
             % Shot Noise
-            shot_noise = sqrt(k * obj.fsimu .* yout)   .* randn(1,size(yout,1));
+            shot_noise = sqrt(k * obj.fsimu .* yout)   .* randn(size(yout,1),1);
 
             yout = yout + shot_noise;
 
@@ -50,7 +50,7 @@ classdef Photodiode
             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));
+            therm_noise = sqrt(therm_noise_pow) .* randn(size(yout,1),1);
 
             yout = yout + therm_noise;
 

Classes/Scope.m

@@ -0,0 +1,193 @@
+classdef Scope
+    %SCOPE Summary of this class goes here
+    %   Detailed explanation goes here
+
+    properties
+        fsimu
+        fadc
+        adcresolution
+        quantbuffer
+
+        rand_samplingdelay %on or off
+        samplingdelay %specifiy a sampling delay
+
+        freq_offset %offset of the sampler
+        samp_jitter %include jitter
+
+        fixed_delay %fix the delay of the filter or use minimal delay for kausal system
+        delay  %specify a fixed delay of the filter
+
+        filtertype
+        lpf_bw
+
+        %during construction
+
+        %during process
+        Nout
+    end
+
+    methods
+        function obj = Scope(options)
+            %SCOPE Construct an instance of this class
+            %   Detailed explanation goes here
+            arguments
+                options.fsimu
+                options.fadc
+                options.adcresolution
+                options.quantbuffer
+
+                options.rand_samplingdelay = 0;
+                options.samplingdelay = 0;
+
+                options.freq_offset = 0;
+                options.samp_jitter = 0;
+
+                options.fixed_delay = 0;
+                options.delay = 0;
+
+                options.filtertype = filtertypes.bessel_bilin;
+                options.lpf_bw = 120e9;
+            end
+
+            obj.fsimu = options.fsimu;
+            obj.fadc = options.fadc;
+            obj.adcresolution = options.adcresolution;
+            obj.quantbuffer = options.quantbuffer;
+
+            obj.rand_samplingdelay = options.rand_samplingdelay; % use a randomized sample delay INSTEAD of samplingdelay
+            obj.samplingdelay = options.samplingdelay; %specifiy a sampling delay
+
+            obj.freq_offset = options.freq_offset; %offset of the sampler
+            obj.samp_jitter = options.samp_jitter; %include jitter in [s]
+
+            obj.fixed_delay = options.fixed_delay; %fix the delay of the filter or use minimal delay for kausal system
+            obj.delay = options.delay;  %specify a fixed delay of the filter
+
+            obj.filtertype = options.filtertype;
+            obj.lpf_bw = options.lpf_bw;
+
+        end
+
+        function yout = process(obj,xin)
+            %METHOD1 Summary of this method goes here
+            %   Detailed explanation goes here
+
+            % apply LPF
+            lpf = Filter('filtdegree',4,"f_cutoff",obj.lpf_bw,"fsamp",obj.fsimu,"filterType",obj.filtertype);
+            yout = lpf.process(xin);
+
+            % sample signal
+            % TODO: implement and test the delays. Also look for delay of
+            % lpf filter
+            yout = obj.sampleSignal(yout);
+
+            % quantize signal
+            yout = obj.quantize(yout);
+
+
+        end
+
+        function jitvec = buildSamplingJitterVector(obj)
+            % Generation of sampling jitter vector for real part of signal
+            jitvec = obj.samp_jitter*randn(obj.Nout,1) ;
+
+            % Avoid exceding the simulation block edges relentlessly!
+            jitvec(1)   = abs(jitvec(1)) ;
+            jitvec(end) = -abs(jitvec(end)) ;
+
+        end
+
+        function sampout = sampleSignal(obj,xin)
+
+            % Initialization sampler and output signal in case of 4real
+
+            Tout = 1 / obj.fadc; % output signal sample period [sec]
+
+            Tsig = length(xin) / obj.fsimu ;      % simulation block period [sec]
+
+            if obj.rand_samplingdelay == 1
+                obj.samplingdelay = Tout*rand ;
+            end
+
+            %get number of delayed samples
+            Ndelay        = floor(obj.samplingdelay*obj.fsimu) ;
+
+            %get fractional sample delay
+            fracdelay     = obj.samplingdelay - Ndelay/obj.fsimu ;
+
+            %calc number of samples after down-sampling no offset
+            obj.Nout  = ceil((Tsig - fracdelay)*obj.fadc) ;
+
+            %build a jitter vector 
+            jitvec = obj.buildSamplingJitterVector();
+
+            % Create the index vector of the sampling instances for the real and imaginary part
+            inx_re = ( (0: Tout : (obj.Nout-1) / obj.fadc)' + jitvec ) * obj.fsimu + 1 ;
+
+            % Create the vector of the relative location of each sample between
+            % to neighboring "Analog" samples for the real and imaginary part
+            val_re = mod(inx_re,1) ;
+
+            % Shift the signal
+            xin = circshift(xin, [0 - Ndelay]) ;
+
+            % Sample the signal using linear interpolation
+            sampout = xin(floor(inx_re)) .* (1 - val_re) + xin(ceil(inx_re)) .* val_re ;
+
+        end
+
+        function quantout = quantize(obj,xin)
+
+            quant_steps = round(2^obj.adcresolution) + rem(round(2^obj.adcresolution),2);
+
+            if quant_steps > 0
+
+                index = round(0.1*length(xin)); %sioe 2022: short blocksize resulted in error -> now 10%-90% of length are evaluated
+
+                max_re = max(real(xin(index:end-index)))*(1+obj.quantbuffer/2); % adaptations in the indices
+                min_re = min(real(xin(index:end-index)))*(1-obj.quantbuffer/2); %(1+para.max/2);
+
+                range_re = max_re-min_re;
+
+                max_im = max(imag(xin))*(1+obj.quantbuffer/2);
+                min_im = min(imag(xin))*(1+obj.quantbuffer/2);
+
+                range_im = max_im-min_im;
+
+                if isreal(xin)
+                    % shift signal and clip to quantizer intervall
+                    xin =   min(max(xin-min_re,0),range_re);
+
+                    % quantize signal
+                    xin = round((quant_steps-1)/range_re*xin);
+
+                    % scale and shift back
+                    quantout = xin*range_re/(quant_steps-1) +min_re;
+
+                else
+
+                    % shift signal and clip to quantizer intervall
+                    xin_re =      min(max(real(xin)-min_re,0),range_re);
+                    xin_im =      min(max(imag(xin)-min_im,0),range_im);
+
+                    % quantize signal
+
+                    % quantize signal and shift back signal
+                    xin_re = round((quant_steps-1)/range_re*xin_re);
+                    xin_im = round((quant_steps-1)/range_im*xin_im);
+
+                    % scale and shift back
+                    quant_out_re = xin_re*range_re/(quant_steps-1) + min_re;
+                    quant_out_im = xin_im*range_im/(quant_steps-1) + min_im;
+                    quantout = quant_out_re + 1i * quant_out_im;
+
+                end
+            else
+                quantout = xin;
+            end
+
+        end
+
+    end
+end
+

Classes/Signal.m

@@ -9,13 +9,14 @@ classdef Signal
     end
     
     methods
-        function obj = Signal(signal, fsym, fsimu)
+        function obj = Signal(signal, fsym)
             %SIGNAL Construct an instance of this class
             %   Detailed explanation goes here
             obj.signal = signal;
             
             obj.fsym = fsym;
-            obj.fsimu = fsimu;
+
+            obj.fsimu = fsym;
 
         end
 
@@ -39,8 +40,8 @@ classdef Signal
             end
         end
 
-        function outputArg = method1(obj,inputArg)
+        function normalizedsignal = signal_(obj)
-
+            normalizedsignal = obj.signal / norm(obj.signal);
         end
 
     end

comm_tb.m

@@ -1,5 +1,5 @@
-O = 10; %order of prbs
+O = 14; %order of prbs
-N = 64; %length of prbs
+N = 2^(O-1); %length of prbs
 [~,seed] = prbs(O,1); %initialize first seed of prbs
 
 % Modulation
@@ -12,6 +12,8 @@ fsym = 56e9;
 fdac = 120e9;
 % Simulation oversampling rate "k";
 kover = 16;
+% ADC Rate
+fadc = 256e9;
 % Simulation frequency in "analog domain"
 fsimu = kover * fdac ;
 
@@ -22,44 +24,73 @@ for i = 1:log2(M)
     [bitpattern(:,i),seed] = prbs(O,N,seed);
 end
 
+S = Signal(bitpattern,fsym);
+
 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);
+fil = Filter('filtdegree',4,"f_cutoff",50e9,"fsamp",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);
+extmodlaser = EML("mode",emlmodes.im_cosinus,"power",10,"fsimu",fsimu,"lambda",1550,"bias",vbias,"u_pi",u_pi,"linewidth",1000000);
-laserfield = eml.process(filtered); 
+laserfield = extmodlaser.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("fsimu",fdac*kover,"fiber_length",20,"alpha",0.2,"D",17,"lambda0",1550);
+fib = Fiber("fsimu",fdac*kover,"fiber_length",0,"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);
 
+fil = Filter('filtdegree',4,"f_cutoff",70e9,"fsamp",fdac,"filterType",filtertypes.bessel_bilin);
+phdiod_out2 = fil.process(phdiod_out);
 
 
+scp = Scope("fsimu",fdac*kover,"fadc",fadc,...
+    "delay",0,"fixed_delay",0,"lpf_bw",120e9,"filtertype",filtertypes.bessel_inp,...
+    "samplingdelay",0,"rand_samplingdelay",0,"freq_offset",0,"samp_jitter",0,...
+    "adcresolution",6,"quantbuffer",0.1);
+
+scope_out = scp.process(phdiod_out2);
+
 figure;
-plot(shapedData);
 hold on
-plot(awgSignal,'DisplayName','skew 0');
+plot( nmlze(resample(shapedData,fadc,fdac)),'DisplayName','shapedData Data');
-plot(filtered,'DisplayName','filtered');
+plot( nmlze(resample(awgSignal,fadc,fdac*kover)),'DisplayName','AWG Data');
-plot(abs(laserfield),'DisplayName','laser');
+plot( nmlze(scope_out),'DisplayName','scope out out');
-plot(abs(att_out),'DisplayName','att_out');
-plot(abs(fib_out),'DisplayName','att_out');
-
 hold off
 
+resample_out = resample(scope_out,fsym,fadc);
+
+
+% figure;
+% hold on
+% 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','fiber_out');
+% plot(abs(phdiod_out),'DisplayName','photo diode out');
+% hold off
+
+% figure;
+% hold on
+% plot(nmlze(pamData),'DisplayName','PAM Data');
+% plot(nmlze(resample_out),'DisplayName','system out');
+% hold off
+
+
+function y = nmlze(x)
+    y = (x-min(x)) / max((x-min(x)));
+end
 
 function pam_sig = pam_mapping(bitpattern, M, unipolar)