diff --git a/Classes/00_signals/Informationsignal.m b/Classes/00_signals/Informationsignal.m index d098d2f..14279f3 100644 --- a/Classes/00_signals/Informationsignal.m +++ b/Classes/00_signals/Informationsignal.m @@ -26,11 +26,11 @@ classdef Informationsignal < Signal end - % function pow = power(obj) - % - % pow = mean(abs(obj.signal.^2),"all") ; - % - % end + function pow = power(obj) + + pow = mean(abs(obj.signal.^2),"all") ; + + end end diff --git a/Classes/00_signals/Signal.m b/Classes/00_signals/Signal.m index 0a45e32..407691c 100644 --- a/Classes/00_signals/Signal.m +++ b/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) diff --git a/Functions/Channel_structures/awgn_alpha_d_channel.m b/Functions/Channel_structures/awgn_alpha_d_channel.m new file mode 100644 index 0000000..4c8375c --- /dev/null +++ b/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 diff --git a/Functions/Channel_structures/awgn_channel.m b/Functions/Channel_structures/awgn_channel.m index 6cf5580..ea03d78 100644 --- a/Functions/Channel_structures/awgn_channel.m +++ b/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; - -end \ No newline at end of file + 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 diff --git a/Functions/EQ_structures/dsp_scope_signal.m b/Functions/EQ_structures/dsp_scope_signal.m index bbb376e..dbe91df 100644 --- a/Functions/EQ_structures/dsp_scope_signal.m +++ b/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 @@ -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 diff --git a/projects/IMDD_base_system/imdd_it.m b/projects/IMDD_base_system/imdd_it.m index 16d715e..ddce315 100644 --- a/projects/IMDD_base_system/imdd_it.m +++ b/projects/IMDD_base_system/imdd_it.m @@ -15,7 +15,7 @@ close all; if 1 - uloops = struct; + uloops = struct; uloops.precomp = [0]; uloops.bitrate = [300].*1e9; %[300,330,360,390,420,450,480] [224,336,360,390,420,448] for MPI % uloops.laser_wavelength = [1293,1297.5,1302,1306.5,1310,1313.4,1318,1322.7,1327.4]; @@ -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.duob_mode = db_mode.db_precoded; + uloops.channel_alpha = [0.5]; + 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 - a=wh.getStoValue('ber',0, [300].*1e9 , 1293, 4, uloops.link_length,alpha,uloops.duob_mode,uloops.decoding_mode); - % ffe = cellfun(@(x) x.ffe_results.metrics.BER, a); - ffe = cellfun(@(x) x.dbt_results.metrics.BER, a); - plot(uloops.link_length,ffe,'DisplayName',sprintf('Alpha: %d',alpha),'LineStyle','-','HandleVisibility','on'); +for alpha = uloops.channel_alpha + 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); + gmi_ffe = cellfun(@(x) x.ffe_results.metrics.GMI, a); + alpha_ffe = cellfun(@(x) x.ffe_results.metrics.Alpha, a); + + 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]); -yline([3.8e-3, 2e-2],'HandleVisibility','off'); +% ylim([5e-5 0.5]); +% yline([3.8e-3, 2e-2],'HandleVisibility','off'); legend beautifyBERplot() ylabel('BER'); diff --git a/projects/IMDD_base_system/imdd_model.m b/projects/IMDD_base_system/imdd_model.m index fdbf366..eb55538 100644 --- a/projects/IMDD_base_system/imdd_model.m +++ b/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,52 +114,74 @@ Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rc","pulselength",1 Digi_sig.spectrum("displayname",'Digi Spectrum','fignum',10,'normalizeTo0dB',1); -%%%%% 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); -% El_sig.spectrum("displayname",'Digi Spectrum','fignum',100,'normalizeTo0dB',0); -% El_sig = El_sig.setPower(0,"dBm"); -%%%%% Low-pass el. components %%%%%% -tx_bwl = tx_bw_nyquist.*f_nyquist; -% tx_bwl = 80e9; -El_sig = Filter('filtdegree',4,"f_cutoff",tx_bwl,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(El_sig); -% El_sig.spectrum("displayname",'Digi Spectrum','fignum',100,'normalizeTo0dB',1); +channel_mode = lower(string(channel_mode)); -%%%%% Electrical Driver Amplifier %%%%%% -% El_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","gain","amplification_db",3).process(El_sig); -El_sig = El_sig.normalize("mode","oneone"); -scaling = 0.6*(u_pi/2-abs(vbias-u_pi/2)); -El_sig = El_sig .* scaling; +switch channel_mode + case "physical" -%%%%% 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); + %%%%% 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); + % El_sig.spectrum("displayname",'Digi Spectrum','fignum',100,'normalizeTo0dB',0); + % El_sig = El_sig.setPower(0,"dBm"); -% Opt_sig.spectrum("displayname",'Opt Spectrum','fignum',10,'normalizeTo0dB',1); + %%%%% Low-pass el. components %%%%%% + tx_bwl = tx_bw_nyquist.*f_nyquist; + % tx_bwl = 80e9; + El_sig = Filter('filtdegree',4,"f_cutoff",tx_bwl,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(El_sig); + % El_sig.spectrum("displayname",'Digi Spectrum','fignum',100,'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); + %%%%% Electrical Driver Amplifier %%%%%% + % El_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","gain","amplification_db",3).process(El_sig); + El_sig = El_sig.normalize("mode","oneone"); + scaling = 0.6*(u_pi/2-abs(vbias-u_pi/2)); + El_sig = El_sig .* scaling; -%%%%%% ROP %%%%%% -Rx_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",rop).process(Opt_sig); + %%%%% 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",eml_alpha).process(El_sig); -%%%%%% PD Square Law %%%%%% -Rx_sig = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20,"nep",1.8e-11).process(Rx_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); -%%%%%% Low-pass RX (PD, El. Connectors and Scope %%%%%% -rx_bwl = rx_bw_nyquist.*f_nyquist; -% rx_bwl = 80e9; -Rx_sig = Filter('filtdegree',4,"f_cutoff",rx_bwl,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(Rx_sig); + %%%%%% ROP %%%%%% + Rx_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",rop).process(Opt_sig); -% %%%%%% Low-pass Scope %%%%%% -Lp_scpe = Filter('filtdegree',4,"f_cutoff",110e9,"fs",fadc,"filterType",filtertypes.butterworth,"active",true); + %%%%%% PD Square Law %%%%%% + Rx_sig = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20,"nep",1.8e-11).process(Rx_sig); -% Rx_sig.spectrum("displayname",'Analog Rx Spectrum','fignum',100,'normalizeTo0dB',1); + %%%%%% Low-pass RX (PD, El. Connectors and Scope %%%%%% + rx_bwl = rx_bw_nyquist.*f_nyquist; + % rx_bwl = 80e9; + Rx_sig = Filter('filtdegree',4,"f_cutoff",rx_bwl,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(Rx_sig); -%%%%%% Scope %%%%%% -Scpe_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); + % %%%%%% Low-pass Scope %%%%%% + Lp_scpe = Filter('filtdegree',4,"f_cutoff",110e9,"fs",fadc,"filterType",filtertypes.butterworth,"active",true); + + % Rx_sig.spectrum("displayname",'Analog Rx Spectrum','fignum',100,'normalizeTo0dB',1); + + %%%%%% Scope %%%%%% + 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 @@ -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, ...