Add AWGN channel models and route IMDD simulation through them

Classes/00_signals/Informationsignal.m

@@ -26,11 +26,11 @@ classdef Informationsignal < Signal
             
         end
 
-        % function pow = power(obj)
+        function pow = power(obj)
-        % 
+
-        %     pow = mean(abs(obj.signal.^2),"all") ;
+            pow = mean(abs(obj.signal.^2),"all") ;
-        % 
+
-        % end
+        end
 
 
     end

Classes/00_signals/Signal.m

@@ -172,9 +172,9 @@ classdef Signal
 
             hold on;
             if isempty(options.color)
-                plot(t* 1e6, sig(1:length(t)), 'DisplayName', dn, 'LineWidth', 0.1, 'Marker', 'none', 'LineStyle','-', 'MarkerSize', 0.1);
+                plot(t* 1e6, sig(1:length(t)), 'DisplayName', dn, 'LineWidth', 0.1, 'Marker', '.', 'LineStyle','none', 'MarkerSize', 0.1);
             else
-                plot(t* 1e6, sig(1:length(t)), 'DisplayName', dn, 'LineWidth', 0.1, 'Marker', 'none', 'LineStyle','-', 'MarkerSize', 0.1,'Color',options.color);
+                plot(t* 1e6, sig(1:length(t)), 'DisplayName', dn, 'LineWidth', 0.1, 'Marker', '.', 'LineStyle','none', 'MarkerSize', 0.1,'Color',options.color);
             end
             % 2 c)
             % - xlabel if not already here: time in readable format (1 ms and not 1e-3 s)

Functions/Channel_structures/awgn_alpha_d_channel.m

@@ -0,0 +1,32 @@
+function signal_out = awgn_alpha_d_channel(signal_in, options)
+%AWGN_ALPHA_D_CHANNEL Apply a 1 + alpha*D FIR channel followed by AWGN.
+
+    arguments
+        signal_in
+        options.alpha (1,1) double = 0
+        options.delay_samples (1,1) double {mustBeInteger, mustBePositive} = 1
+        options.snr_dB (1,1) double = 20
+        options.randkey = []
+    end
+
+    signal_out = signal_in;
+
+    taps = zeros(1, options.delay_samples + 1);
+    taps(1) = 1;
+    taps(end) = options.alpha;
+
+    signal_out.signal = filter(taps, 1, signal_in.signal, [], 1);
+
+    if isa(signal_out, 'Signal')
+        desc = sprintf('1 + alpha*D channel with alpha %.4f and delay %d samples', ...
+            options.alpha, options.delay_samples);
+        meta = struct( ...
+            'alpha', options.alpha, ...
+            'delay_samples', options.delay_samples);
+        signal_out = signal_out.logbookentry(desc, meta);
+    end
+
+    signal_out = awgn_channel(signal_out, ...
+        "snr_dB", options.snr_dB, ...
+        "randkey", options.randkey);
+end

Functions/Channel_structures/awgn_channel.m

@@ -1,5 +1,43 @@
-function signal_out = awgn_channel(signal_in)
+function signal_out = awgn_channel(signal_in, options)
+%AWGN_CHANNEL Apply additive white Gaussian noise to a signal object.
+
+    arguments
+        signal_in
+        options.snr_dB (1,1) double = 20
+        options.randkey = []
+    end
 
     signal_out = signal_in;
+    x = signal_in.signal;
 
+    signalPower = mean(abs(x).^2, 'all'); %or use signal_in.power
+    if signalPower == 0
+        return
+    end
+
+    snrLinear = 10^(options.snr_dB/10);
+    noisePower = signalPower / snrLinear;
+
+    if isempty(options.randkey)
+        if isreal(x)
+            noise = sqrt(noisePower) * randn(size(x));
+        else
+            noise = sqrt(noisePower / 2) * (randn(size(x)) + 1i * randn(size(x)));
+        end
+    else
+        rs = RandStream('mt19937ar', 'Seed', options.randkey);
+        if isreal(x)
+            noise = sqrt(noisePower) * randn(rs, size(x));
+        else
+            noise = sqrt(noisePower / 2) * (randn(rs, size(x)) + 1i * randn(rs, size(x)));
+        end
+    end
+
+    signal_out.signal = x + noise;
+
+    if isa(signal_out, 'Signal')
+        desc = sprintf('AWGN channel with SNR %.2f dB', options.snr_dB);
+        meta = struct('snr_dB', options.snr_dB, 'randkey', options.randkey);
+        signal_out = signal_out.logbookentry(desc, meta);
+    end
 end

Functions/EQ_structures/dsp_scope_signal.m

@@ -47,12 +47,12 @@ function output = dsp_scope_signal(Scpe_sig_raw, Symbols, Tx_bits, options)
     use_dd_mode = 1;
 
     use_ffe         = 1;
-    use_dfe         = 1;
+    use_dfe         = 0;
     use_vnle_mlse   = 0;
     use_dbtgt       = 0;
     use_dbenc       = 0;
     use_ml_mlse     = 0;
-    showAnalysis    = 0;
+    showAnalysis    = 1;
     decoding_mode   = [];
 
     addProcessingResultToDatabase = 0; %#ok<NASGU>
@@ -102,7 +102,7 @@ function output = dsp_scope_signal(Scpe_sig_raw, Symbols, Tx_bits, options)
                 "eth_style_symbol_mapping", 0);
 
             ffe_results.config.equalizer_structure = "ffe";
-            ffe_results.metrics.print;
+            ffe_results.metrics.print("description",'FFE');
             output.ffe_package = ffe_results;
         end
 

projects/IMDD_base_system/imdd_it.m

@@ -23,30 +23,36 @@ if 1
     uloops.M            = [4];
     uloops.link_length  = 1;
     % uloops.link_length  = [0:2:10]; % 1,2,3,5,6,8,10
-    uloops.alpha = [0];
+    uloops.channel_alpha = [0.5];
-    uloops.duob_mode = db_mode.db_precoded;
+    uloops.duob_mode = db_mode.no_db;
     uloops.decoding_mode = "memoryless";
+    uloops.channel_mode = "awgn_alphad";
+    uloops.channel_snr_dB = [20:-2:0];
     
     wh = DataStorage(uloops);
     wh.addStorage("ber");
     
-    wh = submit_handle(@imdd_model,wh,"parallel",0);
+    wh = submit_handle(@imdd_model,wh,"parallel",1);
     
 end
 
 %%
 figure
 hold on
-for alpha = uloops.alpha
+for alpha = uloops.channel_alpha
-    a=wh.getStoValue('ber',0, [300].*1e9 , 1293, 4, uloops.link_length,alpha,uloops.duob_mode,uloops.decoding_mode);
+    a=wh.getStoValue('ber',0, [300].*1e9 , 1293, 4, uloops.link_length,alpha,uloops.duob_mode,uloops.decoding_mode,uloops.channel_mode,uloops.channel_snr_dB);
-    % ffe = cellfun(@(x) x.ffe_results.metrics.BER, a);
+    gmi_ffe = cellfun(@(x) x.ffe_results.metrics.GMI, a);
-    ffe = cellfun(@(x) x.dbt_results.metrics.BER, a);
+    alpha_ffe = cellfun(@(x) x.ffe_results.metrics.Alpha, a);
-    plot(uloops.link_length,ffe,'DisplayName',sprintf('Alpha: %d',alpha),'LineStyle','-','HandleVisibility','on');
+
+    gmi_mlse = cellfun(@(x) x.mlse_results.metrics.GMI, a);
+
+    plot(uloops.channel_snr_dB,gmi_ffe,'DisplayName',sprintf(''),'LineStyle','-','HandleVisibility','on');
+    plot(uloops.channel_snr_dB,gmi_mlse,'DisplayName',sprintf(''),'LineStyle','-','HandleVisibility','on');
 end
 
 set(gca, 'YScale', 'log');
-ylim([5e-5 0.4]);
+% ylim([5e-5 0.5]);
-yline([3.8e-3, 2e-2],'HandleVisibility','off');
+% yline([3.8e-3, 2e-2],'HandleVisibility','off');
 legend
 beautifyBERplot()
 ylabel('BER');

projects/IMDD_base_system/imdd_model.m

@@ -32,6 +32,10 @@ tx_bw_nyquist   = 0.8;
 
 % Channel
 link_length     = 1;
+channel_mode = "physical";
+channel_snr_dB = 20;
+channel_alpha = 0;
+channel_delay_symbols = 1;
 
 % RX
 rop = -8;
@@ -46,7 +50,7 @@ dfe_order = [0 0 0];
 
 pf_ncoeffs = 1;
 
-alpha = 0;
+eml_alpha = 0;
 
 len_tr  = 4096*2;
 
@@ -95,15 +99,12 @@ if fsym_ ~= fsym
 end
 
 
-
-
-
 f_nyquist       = fsym/2;
 
 Pform       = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rc","pulselength",16,"alpha",rcalpha);
 
 [Digi_sig,Symbols,Tx_bits] = PAMsource(...
-    "fsym",fsym,"M",M,"order",18,"useprbs",0,...
+    "fsym",fsym,"M",M,"order",16,"useprbs",0,...
     "fs_out",fdac,...
     "applyclipping",0,"clipfactor",1.5,...
     "applypulseform",apply_pulsef,"pulseformer",Pform,...
@@ -113,6 +114,12 @@ Pform       = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rc","pulselength",1
 
 Digi_sig.spectrum("displayname",'Digi Spectrum','fignum',10,'normalizeTo0dB',1);
 
+
+channel_mode = lower(string(channel_mode));
+
+switch channel_mode
+    case "physical"
+
         %%%%% AWG
         % El_sig = M8199A("kover",kover).process(Digi_sig);
         El_sig = AWG("fdac",fdac,"f_cutoff",fsym,"lpf_active",0,"kover",kover,"bit_resolution",12,"upsampling_method","samplehold","precomp_sinc_rolloff",1).process(Digi_sig);
@@ -132,10 +139,9 @@ scaling = 0.6*(u_pi/2-abs(vbias-u_pi/2));
         El_sig = El_sig .* scaling;
 
         %%%%% MODULATE E/O CONVERSION %%%%%%
-[Opt_sig] = EML("mode",eml_mode.im_cosinus,"power",3,"fsimu",El_sig.fs,"lambda",laser_wavelength,"bias",vbias,"u_pi",u_pi,"linewidth",laser_linewidth,"randomkey",random_key+1,"alpha",alpha).process(El_sig);
+        [Opt_sig] = EML("mode",eml_mode.im_cosinus,"power",3,"fsimu",El_sig.fs,"lambda",laser_wavelength,"bias",vbias,"u_pi",u_pi,"linewidth",laser_linewidth,"randomkey",random_key+1,"alpha",eml_alpha).process(El_sig);
 
         % Opt_sig.spectrum("displayname",'Opt Spectrum','fignum',10,'normalizeTo0dB',1);
-
         Opt_sig    = Fiber("fsimu",Opt_sig.fs,"fiber_length",link_length,"alpha",0.3,"D",0,"lambda0",1310,"gamma",0,"Dslope",0.07).process(Opt_sig);
 
         %%%%%% ROP %%%%%%
@@ -155,11 +161,28 @@ Lp_scpe = Filter('filtdegree',4,"f_cutoff",110e9,"fs",fadc,"filterType",filterty
         % Rx_sig.spectrum("displayname",'Analog Rx Spectrum','fignum',100,'normalizeTo0dB',1);
 
         %%%%%% Scope %%%%%%
-Scpe_sig = Scope("fsimu",fdac*kover,"fadc",fadc,...
+        Rx_sig = Scope("fsimu",fdac*kover,"fadc",fadc,...
             "delay",0,"fixed_delay",0,"filtertype",filtertypes.butterworth,...
             "samplingdelay",0,"rand_samplingdelay",0,"freq_offset",0,"samp_jitter",0,...
             "adcresolution",8,"quantbuffer",0.1,'block_dc',1,'lpf_active',1,'H_lpf',Lp_scpe).process(Rx_sig);
 
+    case "awgn"
+        Rx_sig = awgn_channel(Digi_sig, ...
+            "snr_dB", channel_snr_dB, ...
+            "randkey", random_key + 2);
+
+    case {"awgn_alphad", "awgn_1alphad", "awgn_1_plus_alpha_d"}
+        delaySamples = max(1, round(channel_delay_symbols * Digi_sig.fs / fsym));
+        Rx_sig = awgn_alpha_d_channel(Digi_sig, ...
+            "alpha", channel_alpha, ...
+            "delay_samples", delaySamples, ...
+            "snr_dB", channel_snr_dB, ...
+            "randkey", random_key + 2);
+
+    otherwise
+        error('Unknown channel_mode "%s". Supported: physical, awgn, awgn_alphaD.', channel_mode);
+end
+
 txPulseformer = [];
 if apply_pulsef
     txPulseformer = Pform;
@@ -171,9 +194,9 @@ dspParameters.mu_ffe = mu_ffe;
 dspParameters.mu_dfe = mu_dfe;
 dspParameters.mu_dc = mu_dc;
 dspParameters.use_ffe = 1;
-dspParameters.use_dfe = 1;
+dspParameters.use_dfe = 0;
 dspParameters.use_vnle_mlse = 1;
-dspParameters.use_dbtgt = 1;
+dspParameters.use_dbtgt = 0;
 dspParameters.use_dbenc = 0;
 dspParameters.use_ml_mlse = 0;
 dspParameters.showAnalysis = 0;
@@ -187,7 +210,7 @@ dspParameters.ffe_order_dbtgt = [50, 0, 0];
 dspParameters.dfe_order_dbtgt = dfe_order;
 dspParameters.decoding_mode = decoding_mode;
 
-scopeOutput = dsp_scope_signal(Scpe_sig, Symbols, Tx_bits, ...
+scopeOutput = dsp_scope_signal(Rx_sig, Symbols, Tx_bits, ...
     "fsym", fsym, ...
     "M", M, ...
     "duob_mode", duob_mode, ...