Merge branch 'main' of cau-git.rz.uni-kiel.de:nt/mitarbeiter/silas/imdd_simulation
This commit is contained in:
Silas Oettinghaus
@@ -1,153 +1,187 @@
|
||||
function [db_results] = duobinary_target(eq_, mlse_,M, rx_signal, tx_symbols, tx_bits, options)
|
||||
|
||||
arguments
|
||||
eq_
|
||||
mlse_
|
||||
M
|
||||
rx_signal
|
||||
tx_symbols
|
||||
tx_bits
|
||||
options.precode_mode db_mode
|
||||
options.showAnalysis = 0;
|
||||
options.eth_style_symbol_mapping = 0;
|
||||
options.postFFE = [];
|
||||
end
|
||||
arguments
|
||||
eq_
|
||||
mlse_
|
||||
M
|
||||
rx_signal
|
||||
tx_symbols
|
||||
tx_bits
|
||||
options.precode_mode db_mode
|
||||
options.showAnalysis = 0;
|
||||
options.eth_style_symbol_mapping = 0;
|
||||
options.postFFE = [];
|
||||
options.decoding_mode db_decoder = db_decoder.sequencedetection;
|
||||
end
|
||||
|
||||
|
||||
%Duobinary Targeting
|
||||
db_ref_sequence = Duobinary().encode(tx_symbols);
|
||||
db_ref_constellation = unique(db_ref_sequence.signal);
|
||||
[eq_signal, eq_noise] = eq_.process(rx_signal,db_ref_sequence);
|
||||
|
||||
if ~isempty(options.postFFE)
|
||||
[eq_signal,eq_noise] = options.postFFE.process(eq_signal,db_ref_sequence);
|
||||
end
|
||||
%
|
||||
|
||||
switch options.decoding_mode
|
||||
case db_decoder.sequencedetection %MLSE
|
||||
mlse_.DIR = [1,1];
|
||||
if isa(mlse_,'MLSE_viterbi')
|
||||
mlse_sig_sd = mlse_.process(eq_signal);
|
||||
else
|
||||
[mlse_sig_sd,LLR,GMI_MLSE] = mlse_.process(eq_signal,tx_symbols);
|
||||
end
|
||||
pam_sig_hd = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).quantize(mlse_sig_sd);
|
||||
case db_decoder.memoryless %DB Target FFE
|
||||
% Hard decision on FFE output
|
||||
eq_signal_hd = PAMmapper(M, 0).quantize(eq_signal,'custom_const',db_ref_constellation.');
|
||||
eq_signal_hd = Duobinary().decode(eq_signal_hd);
|
||||
pam_sig_hd = eq_signal_hd;
|
||||
% tx_symbols = Duobinary().encode(tx_symbols);
|
||||
% tx_symbols = Duobinary().decode(tx_symbols.*db_const_scale_factor);
|
||||
end
|
||||
|
||||
% precoding to mitigate error propagation, most prominently used in
|
||||
% combination with duobinary signaling to avoid catastrophic error
|
||||
% behavior (see J.W.M. Bergmans, Digital Baseband Transmission and Recording -> partial response signaling)
|
||||
switch options.precode_mode
|
||||
|
||||
case db_mode.no_db
|
||||
% TX Data is not precoded:
|
||||
|
||||
% A) Emulate diff precoding
|
||||
switch options.decoding_mode
|
||||
case db_decoder.sequencedetection %MLSE
|
||||
mlse_sig_hd_precoded = Duobinary().encode(pam_sig_hd);
|
||||
mlse_sig_hd_precoded = Duobinary().decode(mlse_sig_hd_precoded);
|
||||
case db_decoder.memoryless %DB Target FFE
|
||||
% mlse_sig_hd_precoded = Duobinary().encode(pam_sig_hd);
|
||||
% mlse_sig_hd_precoded = Duobinary().decode(mlse_sig_hd_precoded);
|
||||
mlse_sig_hd_precoded = pam_sig_hd;
|
||||
end
|
||||
|
||||
%Duobinary Targeting
|
||||
db_ref_sequence = Duobinary().encode(tx_symbols);
|
||||
db_ref_constellation = unique(db_ref_sequence.signal);
|
||||
[eq_signal, eq_noise] = eq_.process(rx_signal,db_ref_sequence);
|
||||
tx_symbols_precoded = Duobinary().encode(tx_symbols);
|
||||
tx_symbols_precoded = Duobinary().decode(tx_symbols_precoded);
|
||||
|
||||
if ~isempty(options.postFFE)
|
||||
[eq_signal,eq_noise] = options.postFFE.process(eq_signal,db_ref_sequence);
|
||||
end
|
||||
tx_bits_precoded = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(tx_symbols_precoded);
|
||||
rx_bits_mlse = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(mlse_sig_hd_precoded);
|
||||
|
||||
[~,errors_db_diff_precoded,ber_db_diff_precoded,~] = calc_ber(rx_bits_mlse.signal,tx_bits_precoded.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
|
||||
|
||||
%B) Just determine BER
|
||||
rx_bits_mlse = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(pam_sig_hd);
|
||||
tx_bits = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(tx_symbols);
|
||||
[bits_mlse,errors_mlse,ber_db,~] = calc_ber(rx_bits_mlse.signal,tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
|
||||
|
||||
case db_mode.db_precoded
|
||||
|
||||
% Daten SIND TATSÄCHLICH precoded auf TX Seite:
|
||||
|
||||
% A) Decode at Rx if no DB targeting was applied (we are in VNLE or MLSE EQ structure here!
|
||||
switch options.decoding_mode
|
||||
case db_decoder.sequencedetection %MLSE
|
||||
mlse_sig_hd_decoded = Duobinary().encode(pam_sig_hd,"M",M);
|
||||
mlse_sig_hd_decoded = Duobinary().decode(mlse_sig_hd_decoded,"M",M);
|
||||
case db_decoder.memoryless %DB Target FFE
|
||||
mlse_sig_hd_decoded = pam_sig_hd;
|
||||
end
|
||||
|
||||
mlse_.DIR = [1,1];
|
||||
%
|
||||
tx_symbols_precoded = Duobinary().encode(tx_symbols,"M",M);
|
||||
tx_symbols_precoded = Duobinary().decode(tx_symbols_precoded,"M",M);
|
||||
|
||||
if isa(mlse_,'MLSE_viterbi')
|
||||
mlse_sig_sd = mlse_.process(eq_signal);
|
||||
else
|
||||
[mlse_sig_sd,LLR,GMI_MLSE] = mlse_.process(eq_signal,tx_symbols);
|
||||
end
|
||||
tx_bits_precoded = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(tx_symbols_precoded);
|
||||
rx_bits_mlse_decoded = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(mlse_sig_hd_decoded);
|
||||
|
||||
mlse_sig_hd = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).quantize(mlse_sig_sd);
|
||||
|
||||
% precoding to mitigate error propagation, most prominently used in
|
||||
% combination with duobinary signaling to avoid catastrophic error
|
||||
% behavior (see J.W.M. Bergmans, Digital Baseband Transmission and Recording -> partial response signaling)
|
||||
switch options.precode_mode
|
||||
|
||||
case db_mode.no_db
|
||||
% TX Data is not precoded:
|
||||
|
||||
% A) Emulate diff precoding
|
||||
mlse_sig_hd_precoded = Duobinary().encode(mlse_sig_hd,"M",M);
|
||||
mlse_sig_hd_precoded = Duobinary().decode(mlse_sig_hd_precoded,"M",M);
|
||||
|
||||
tx_symbols_precoded = Duobinary().encode(tx_symbols);
|
||||
tx_symbols_precoded = Duobinary().decode(tx_symbols_precoded);
|
||||
|
||||
tx_bits_precoded = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(tx_symbols_precoded);
|
||||
rx_bits_mlse = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(mlse_sig_hd_precoded);
|
||||
|
||||
[~,errors_db_diff_precoded,ber_db_diff_precoded,~] = calc_ber(rx_bits_mlse.signal,tx_bits_precoded.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
|
||||
|
||||
%B) Just determine BER
|
||||
rx_bits_mlse = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(mlse_sig_hd);
|
||||
[bits_mlse,errors_mlse,ber_db,~] = calc_ber(rx_bits_mlse.signal,tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
|
||||
|
||||
case db_mode.db_precoded
|
||||
|
||||
% Daten SIND TATSÄCHLICH precoded auf TX Seite:
|
||||
|
||||
% A) Decode at Rx if no DB targeting was applied (we are in VNLE or MLSE EQ structure here!
|
||||
mlse_sig_hd_decoded = Duobinary().encode(mlse_sig_hd,"M",M);
|
||||
mlse_sig_hd_decoded = Duobinary().decode(mlse_sig_hd_decoded,"M",M);
|
||||
rx_bits_mlse_decoded = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(mlse_sig_hd_decoded);
|
||||
[~,errors_db_diff_precoded,ber_db_diff_precoded,a] = calc_ber(rx_bits_mlse_decoded.signal,tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
|
||||
burst_db_precoded = count_error_bursts(a, 40);
|
||||
% B) Omit the Coding by comparing with demapped TX symbol sequence
|
||||
|
||||
tx_bits = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(tx_symbols);
|
||||
rx_bits_mlse = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(mlse_sig_hd);
|
||||
[bits_db,errors_db,ber_db,a] = calc_ber(rx_bits_mlse.signal,tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
|
||||
burst_db = count_error_bursts(a, 40);
|
||||
|
||||
cols = linspecer(8);
|
||||
figure();hold on;
|
||||
stem(1:40,burst_db,'LineWidth',1,'Color',cols(4,:),'Marker','_','DisplayName','w/o diff. precoder');
|
||||
stem(1:40,burst_db_precoded,'LineWidth',1,'Color',cols(3,:),'Marker','.','LineStyle','-','DisplayName','w diff. precoder');
|
||||
xlabel('Bit Error Burst Length')
|
||||
ylabel('Occurence')
|
||||
set(gca, 'yscale', 'log');
|
||||
end
|
||||
|
||||
% M = numel(unique(tx_symbols.signal));
|
||||
rx_bits = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(mlse_sig_hd);
|
||||
|
||||
[bits_db,errors_db,ber_db,errorIndice_db] = calc_ber(rx_bits.signal,tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
alpha = arburg(eq_noise.signal,1);%pf_.coefficients(2);
|
||||
alpha = alpha(2);
|
||||
if isa(mlse_,'MLSE_viterbi')
|
||||
gmi_mlse = NaN;
|
||||
air_mlse = NaN;
|
||||
else
|
||||
gmi_mlse = GMI_MLSE;
|
||||
air_mlse = tx_symbols.fs .* floor(log2(double(M))*10)/10 .* gmi_mlse ./ log2(double(M));
|
||||
end
|
||||
|
||||
db_results = struct();
|
||||
db_results.metrics = Metricstruct;
|
||||
db_results.metrics.result_id = NaN;
|
||||
db_results.metrics.run_id = NaN;
|
||||
db_results.metrics.eqParam_id = NaN;
|
||||
db_results.metrics.date_of_processing = datetime('now');
|
||||
db_results.metrics.BER = ber_db;
|
||||
db_results.metrics.numBits = bits_db;
|
||||
db_results.metrics.numBitErr = errors_db;
|
||||
db_results.metrics.BER_precoded = ber_db_diff_precoded;
|
||||
db_results.metrics.numBitErr_precoded = errors_db_diff_precoded;
|
||||
db_results.metrics.GMI = gmi_mlse;
|
||||
db_results.metrics.AIR = air_mlse;
|
||||
db_results.metrics.MLSE_dir = mlse_.DIR;
|
||||
db_results.metrics.Alpha = alpha;
|
||||
|
||||
% Create DB results structure
|
||||
db_results.config = Equalizerstruct();
|
||||
eq_.e = [];
|
||||
eq_.e2 = [];
|
||||
eq_.e3 = [];
|
||||
db_results.config.eq = jsonencode(eq_);
|
||||
% mlse_.DIR = [];
|
||||
db_results.config.mlse = jsonencode(mlse_);
|
||||
db_results.config.equalizer_structure = int32(equalizer_structure.vnle_db_mlse);
|
||||
db_results.config.comment = 'function: Duobinary tgt. (VNLE -> MLSE)';
|
||||
|
||||
if options.showAnalysis
|
||||
|
||||
eq_signal.eye(eq_signal.fs,M,"fignum",249);
|
||||
|
||||
|
||||
eq_noise = eq_noise - mean(eq_noise.signal);
|
||||
|
||||
rx_signal.spectrum("normalizeTo0dB",1,"fignum",250,"displayname","Rx Spectrum");
|
||||
|
||||
Duobinary().encode(tx_symbols).spectrum("normalizeTo0dB",1,"fignum",10,"displayname","DB encoded reference");
|
||||
|
||||
showEQNoisePSD(eq_noise,"fignum",250,"displayname",'Duobinary Target Noise after Equalization');
|
||||
|
||||
fprintf('DB tgt BER: %.2e \n',ber_db);
|
||||
|
||||
figure(341); clf;
|
||||
showLevelHistogram(eq_signal, db_ref_sequence, "fignum", 341);
|
||||
|
||||
end
|
||||
[~,errors_db_diff_precoded,ber_db_diff_precoded,a] = calc_ber(rx_bits_mlse_decoded.signal,tx_bits_precoded.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
|
||||
burst_db_precoded = count_error_bursts(a, 40);
|
||||
% B) Omit the Coding by comparing with demapped TX symbol sequence
|
||||
|
||||
tx_bits = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(tx_symbols);
|
||||
rx_bits_mlse = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(pam_sig_hd);
|
||||
[bits_db,errors_db,ber_db,a] = calc_ber(rx_bits_mlse.signal,tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
|
||||
burst_db = count_error_bursts(a, 40);
|
||||
|
||||
cols = linspecer(8);
|
||||
figure();hold on;
|
||||
stem(1:40,burst_db,'LineWidth',1,'Color',cols(4,:),'Marker','_','DisplayName','w/o diff. precoder');
|
||||
stem(1:40,burst_db_precoded,'LineWidth',1,'Color',cols(3,:),'Marker','.','LineStyle','-','DisplayName','w diff. precoder');
|
||||
xlabel('Bit Error Burst Length')
|
||||
ylabel('Occurence')
|
||||
set(gca, 'yscale', 'log');
|
||||
end
|
||||
|
||||
% M = numel(unique(tx_symbols.signal));
|
||||
rx_bits = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(pam_sig_hd);
|
||||
|
||||
[bits_db,errors_db,ber_db,errorIndice_db] = calc_ber(rx_bits.signal,tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
|
||||
|
||||
alpha = arburg(eq_noise.signal,1);%pf_.coefficients(2);
|
||||
alpha = alpha(2);
|
||||
|
||||
switch options.decoding_mode
|
||||
case db_decoder.sequencedetection %MLSE
|
||||
if isa(mlse_,'MLSE_viterbi')
|
||||
gmi_mlse = NaN;
|
||||
air_mlse = NaN;
|
||||
else
|
||||
gmi_mlse = GMI_MLSE;
|
||||
air_mlse = tx_symbols.fs .* floor(log2(double(M))*10)/10 .* gmi_mlse ./ log2(double(M));
|
||||
end
|
||||
case db_decoder.memoryless %DB Target FFE
|
||||
% [gmi] = calc_air(eq_signal_sd, tx_symbols, "skip_front", 10000, "skip_end", 10000);
|
||||
[gmi] = calc_ngmi(eq_signal,tx_symbols);
|
||||
gmi_mlse = max(gmi,0);
|
||||
air_mlse = tx_symbols.fs .* floor(log2(double(M))*10)/10 .* gmi ./ log2(double(M));
|
||||
end
|
||||
|
||||
db_results = struct();
|
||||
db_results.metrics = Metricstruct;
|
||||
db_results.metrics.result_id = NaN;
|
||||
db_results.metrics.run_id = NaN;
|
||||
db_results.metrics.eqParam_id = NaN;
|
||||
db_results.metrics.date_of_processing = datetime('now');
|
||||
db_results.metrics.BER = ber_db;
|
||||
db_results.metrics.numBits = bits_db;
|
||||
db_results.metrics.numBitErr = errors_db;
|
||||
db_results.metrics.BER_precoded = ber_db_diff_precoded;
|
||||
db_results.metrics.numBitErr_precoded = errors_db_diff_precoded;
|
||||
db_results.metrics.GMI = gmi_mlse;
|
||||
db_results.metrics.AIR = air_mlse;
|
||||
db_results.metrics.MLSE_dir = mlse_.DIR;
|
||||
db_results.metrics.Alpha = alpha;
|
||||
|
||||
% Create DB results structure
|
||||
db_results.config = Equalizerstruct();
|
||||
eq_.e = [];
|
||||
eq_.e2 = [];
|
||||
eq_.e3 = [];
|
||||
db_results.config.eq = jsonencode(eq_);
|
||||
% mlse_.DIR = [];
|
||||
db_results.config.mlse = jsonencode(mlse_);
|
||||
db_results.config.equalizer_structure = int32(equalizer_structure.vnle_db_mlse);
|
||||
db_results.config.comment = 'function: Duobinary tgt. (VNLE -> MLSE)';
|
||||
|
||||
if options.showAnalysis
|
||||
|
||||
eq_signal.eye(eq_signal.fs,M,"fignum",249);
|
||||
|
||||
|
||||
eq_noise = eq_noise - mean(eq_noise.signal);
|
||||
|
||||
rx_signal.spectrum("normalizeTo0dB",1,"fignum",250,"displayname","Rx Spectrum");
|
||||
|
||||
Duobinary().encode(tx_symbols).spectrum("normalizeTo0dB",1,"fignum",10,"displayname","DB encoded reference");
|
||||
|
||||
showEQNoisePSD(eq_noise,"fignum",250,"displayname",'Duobinary Target Noise after Equalization');
|
||||
|
||||
fprintf('DB tgt BER: %.2e \n',ber_db);
|
||||
|
||||
figure(341); clf;
|
||||
showLevelHistogram(eq_signal, db_ref_sequence, "fignum", 341);
|
||||
|
||||
end
|
||||
|
||||
|
||||
end
|
||||
211
Functions/EQ_structures/ffe_db.m
Normal file
211
Functions/EQ_structures/ffe_db.m
Normal file
@@ -0,0 +1,211 @@
|
||||
function [ffe_results] = ffe_db(eq_, M, rx_signal, tx_symbols, tx_bits, options)
|
||||
% FFE Processes signals through FFE equalizer
|
||||
%
|
||||
% Inputs:
|
||||
% eq_ - Equalizer object
|
||||
% M - Modulation order
|
||||
% rx_signal - Received signal
|
||||
% tx_symbols - Transmitted symbols
|
||||
% tx_bits - Transmitted bits
|
||||
% options - Optional parameters
|
||||
%
|
||||
% Outputs:
|
||||
% ffe_results - Results from FFE processing
|
||||
|
||||
arguments
|
||||
eq_
|
||||
M
|
||||
rx_signal
|
||||
tx_symbols
|
||||
tx_bits
|
||||
options.db_target = 0;
|
||||
options.precode_mode db_mode
|
||||
options.showAnalysis = 0;
|
||||
options.eth_style_symbol_mapping = 0;
|
||||
options.postFFE = [];
|
||||
options.database = [];
|
||||
end
|
||||
|
||||
%% Process signals through equalizer
|
||||
% FFE or VNLE
|
||||
if options.db_target
|
||||
tx_symbols_ref = Duobinary().encode(tx_symbols);
|
||||
db_ref_constellation = unique(tx_symbols_ref.signal);
|
||||
[eq_signal_sd, eq_noise] = eq_.process(rx_signal, tx_symbols_ref);
|
||||
else
|
||||
[eq_signal_sd, eq_noise] = eq_.process(rx_signal, tx_symbols);
|
||||
end
|
||||
|
||||
% Apply post-FFE if provided
|
||||
if ~isempty(options.postFFE)
|
||||
tic
|
||||
[eq_signal_sd, eq_noise] = options.postFFE.process(eq_signal_sd, tx_symbols);
|
||||
toc
|
||||
end
|
||||
|
||||
try
|
||||
ch_coefficients = arburg(eq_noise.signal,1);
|
||||
channel_alpha = ch_coefficients(2);
|
||||
end
|
||||
|
||||
% Hard decision on FFE output
|
||||
if options.db_target
|
||||
eq_signal_hd = PAMmapper(M, 0).quantize(eq_signal_sd,'custom_const',db_ref_constellation.');
|
||||
else
|
||||
eq_signal_hd = PAMmapper(M, 0).quantize(eq_signal_sd);
|
||||
end
|
||||
|
||||
if options.db_target
|
||||
eq_signal_hd = Duobinary().decode(eq_signal_hd);
|
||||
tx_symbols = Duobinary().encode(tx_symbols);
|
||||
tx_symbols = Duobinary().decode(tx_symbols);
|
||||
end
|
||||
|
||||
%% Calculate BER based on precoding mode
|
||||
[bits, errors, ber, error_pos, errors_precoded, ber_precoded] = calculateBER(eq_signal_hd, tx_symbols, tx_bits, options.precode_mode, M, options.eth_style_symbol_mapping);
|
||||
|
||||
%% Calculate performance metrics
|
||||
[snr, snr_lvl] = calc_snr(tx_symbols.signal, eq_noise.signal);
|
||||
% [gmi] = calc_air(eq_signal_sd, tx_symbols, "skip_front", 10000, "skip_end", 10000);
|
||||
[gmi] = calc_ngmi(eq_signal_sd,tx_symbols);
|
||||
gmi = max(gmi,0);
|
||||
|
||||
air = tx_symbols.fs .* floor(log2(double(M))*10)/10 .* gmi ./ log2(double(M));
|
||||
[evm_total, evm_lvl] = calc_evm(eq_signal_sd, tx_symbols);
|
||||
[std_total, std_lvl] = calc_std(eq_signal_sd, tx_symbols);
|
||||
[std_rxraw_total, std_rxraw_lvl] = calc_std(rx_signal.resample("fs_out", tx_symbols.fs), tx_symbols);
|
||||
|
||||
%% Display analysis if requested
|
||||
if options.showAnalysis
|
||||
displayAnalysis(eq_noise, eq_signal_sd, rx_signal, eq_, tx_symbols, M, options.postFFE);
|
||||
end
|
||||
|
||||
|
||||
%% Prepare output structure
|
||||
% Determine postFFE order
|
||||
if ~isempty(options.postFFE)
|
||||
npostFFE = options.postFFE.order;
|
||||
else
|
||||
npostFFE = 0;
|
||||
end
|
||||
|
||||
% Create FFE results structure
|
||||
ffe_results = struct();
|
||||
try
|
||||
eq_.e = [];
|
||||
eq_.e2 = [];
|
||||
eq_.e3 = [];
|
||||
eq_.b = [];
|
||||
eq_.b2 = [];
|
||||
eq_.b3 = [];
|
||||
end
|
||||
|
||||
ffe_results.config = Equalizerstruct();
|
||||
ffe_results.config.eq = jsonencode(eq_);
|
||||
ffe_results.config.equalizer_structure = int32(equalizer_structure.ffe);
|
||||
ffe_results.config.comment = 'function: ffe';
|
||||
|
||||
ffe_results.metrics = Metricstruct;
|
||||
ffe_results.metrics.result_id = NaN;
|
||||
ffe_results.metrics.run_id = NaN;
|
||||
ffe_results.metrics.eqParam_id = NaN;
|
||||
ffe_results.metrics.date_of_processing = datetime('now');
|
||||
ffe_results.metrics.BER = ber;
|
||||
ffe_results.metrics.numBits = bits;
|
||||
ffe_results.metrics.numBitErr = errors;
|
||||
ffe_results.metrics.BER_precoded = ber_precoded;
|
||||
ffe_results.metrics.numBitErr_precoded = errors_precoded;
|
||||
ffe_results.metrics.SNR = snr;
|
||||
ffe_results.metrics.SNR_level = snr_lvl;
|
||||
ffe_results.metrics.STD = std_total;
|
||||
ffe_results.metrics.STD_level = std_lvl;
|
||||
ffe_results.metrics.STDrx = std_rxraw_total;
|
||||
ffe_results.metrics.STDrx_level = std_rxraw_lvl;
|
||||
ffe_results.metrics.GMI = gmi;
|
||||
ffe_results.metrics.AIR = air;
|
||||
ffe_results.metrics.EVM = evm_total;
|
||||
ffe_results.metrics.EVM_level = evm_lvl;
|
||||
ffe_results.metrics.Alpha = channel_alpha;
|
||||
|
||||
|
||||
end
|
||||
|
||||
%% Helper Functions
|
||||
function [bits, errors, ber, error_pos, errors_precoded, ber_precoded] = calculateBER(eq_signal_hd, tx_symbols, tx_bits, precode_mode, M, eth_style)
|
||||
% Calculate BER based on precoding mode
|
||||
mapper = PAMmapper(M, 0, "eth_style", eth_style);
|
||||
|
||||
switch precode_mode
|
||||
case db_mode.no_db
|
||||
% TX Data is not precoded
|
||||
% A) Emulate diff precoding
|
||||
eq_signal_hd_precoded = Duobinary().encode(eq_signal_hd, "M", M);
|
||||
eq_signal_hd_precoded = Duobinary().decode(eq_signal_hd_precoded, "M", M);
|
||||
|
||||
tx_symbols_precoded = Duobinary().encode(tx_symbols);
|
||||
tx_symbols_precoded = Duobinary().decode(tx_symbols_precoded);
|
||||
|
||||
tx_bits_precoded = mapper.demap(tx_symbols_precoded);
|
||||
|
||||
rx_bits = mapper.demap(eq_signal_hd_precoded);
|
||||
[~, errors_precoded, ber_precoded, ~] = calc_ber(rx_bits.signal, tx_bits_precoded.signal, "skip_front", 10, "skip_end", 10, "returnErrorLocation", 1);
|
||||
|
||||
% B) Just determine BER
|
||||
rx_bits = mapper.demap(eq_signal_hd);
|
||||
tx_bits = mapper.demap(tx_symbols);
|
||||
[bits, errors, ber, error_pos] = calc_ber(rx_bits.signal, tx_bits.signal, "skip_front", 10, "skip_end", 10, "returnErrorLocation", 1);
|
||||
|
||||
case db_mode.db_precoded
|
||||
% Data is precoded on TX side
|
||||
% A) Decode at Rx if no DB targeting was applied
|
||||
eq_signal_hd_decoded = Duobinary().encode(eq_signal_hd, "M", M);
|
||||
eq_signal_hd_decoded = Duobinary().decode(eq_signal_hd_decoded, "M", M);
|
||||
rx_bits_decoded = mapper.demap(eq_signal_hd_decoded);
|
||||
[~, errors_precoded, ber_precoded, ~] = calc_ber(rx_bits_decoded.signal, tx_bits.signal, "skip_front", 10, "skip_end", 10, "returnErrorLocation", 1);
|
||||
|
||||
% B) Omit the Coding by comparing with demapped TX symbol sequence
|
||||
tx_bits_demapped = mapper.demap(tx_symbols);
|
||||
rx_bits = mapper.demap(eq_signal_hd);
|
||||
[bits, errors, ber, error_pos] = calc_ber(rx_bits.signal, tx_bits_demapped.signal, "skip_front", 10, "skip_end", 10, "returnErrorLocation", 1);
|
||||
end
|
||||
end
|
||||
|
||||
function displayAnalysis(eq_noise, eq_signal_sd, rx_signal, eq_, tx_symbols, M, postFFE)
|
||||
% Display analysis plots and metrics
|
||||
|
||||
% Initialize figure handles
|
||||
% Corrected line - added tx_symbols as second positional argument
|
||||
% showLevelScatter(rx_signal.resample("fs_out", tx_symbols.fs), tx_symbols, "fignum", 100);
|
||||
|
||||
warning off
|
||||
showLevelScatter(eq_signal_sd, tx_symbols, "fignum", 101);
|
||||
figure(gcf);hold on; plot(((1:length(eq_noise.signal)) / eq_noise.fs) * 1e6,movmean(eq_noise.signal,2000,1), 'LineWidth',3,'Color','black')
|
||||
warning on
|
||||
|
||||
showLevelHistogram(eq_signal_sd, tx_symbols, "fignum", 102);
|
||||
|
||||
showEQNoisePSD(eq_noise, "fignum", 103, "displayname", 'Residual Noise after FFE');
|
||||
|
||||
% Figure 2: Post-FFE coefficients (if available)
|
||||
if ~isempty(postFFE)
|
||||
showEQcoefficients('n1', postFFE.e, "displayname", 'Coefficients', 'fignum', 104);
|
||||
end
|
||||
|
||||
try
|
||||
figure(339);hold on
|
||||
showEQfilter(eq_.e, eq_signal_sd.fs.*2,"displayname",'training','fignum',339);
|
||||
% showEQfilter(eq_.e, eq_signal_sd.fs.*2,"displayname",'dec. directed','fignum',339);
|
||||
legend on
|
||||
end
|
||||
|
||||
% try
|
||||
% figure(240); hold on; plot(pow2db(movmean(eq_.debug_struct.error_tr',100)));ylim([-30,3]);title('error training');
|
||||
%
|
||||
% figure(241); hold on; plot(pow2db(movmean(eq_.debug_struct.update_tr',100)));title('update step training');
|
||||
%
|
||||
% figure(242); hold on; plot(pow2db(movmean(eq_.debug_struct.update',1000)));title('update step dd');
|
||||
% end
|
||||
|
||||
eq_signal_sd.eye(eq_signal_sd.fs,M,"displayname",'Eye','fignum',105);
|
||||
|
||||
end
|
||||
@@ -9,31 +9,35 @@
|
||||
% S0 - dispersion slope at ZDW [ps/(nm^2·km)]
|
||||
% ------------------------------------------------------------
|
||||
|
||||
lambda0_nm = 1310;
|
||||
S0 = [0.07,0.09]';
|
||||
%% ZDW 1
|
||||
lambda0_nm = [1310];
|
||||
S0 = [0.075,0.095]';
|
||||
% wavelength axis around ZDW
|
||||
lambda_nm = linspace(lambda0_nm-60, lambda0_nm+60, 400);
|
||||
lambda_nm = linspace(1240, 1360, 400);
|
||||
|
||||
% exact dispersion model
|
||||
D_exact = (S0./4) .* ( ...
|
||||
D_exact_1 = (S0./4) .* ( ...
|
||||
lambda_nm - (lambda0_nm^4)./(lambda_nm.^3) );
|
||||
|
||||
% linearized model around ZDW
|
||||
D_lin = S0 .* (lambda_nm - lambda0_nm);
|
||||
%% ZDW 2
|
||||
lambda0_nm = [1320];
|
||||
% exact dispersion model
|
||||
D_exact_2 = (S0./4) .* ( ...
|
||||
lambda_nm - (lambda0_nm^4)./(lambda_nm.^3) );
|
||||
|
||||
% plot
|
||||
|
||||
%% plot
|
||||
figure('Color','w'); hold on;
|
||||
cols = [0.6510 0.8078 0.8902;...
|
||||
0.1216 0.4706 0.7059;...
|
||||
0.6980 0.8745 0.5412;...
|
||||
0.2000 0.6275 0.1725];
|
||||
|
||||
% plot(lambda_nm, D_lin(1,:), '-', 'LineWidth',2, 'DisplayName','Linearized model','Color',[cols(1,:)]);
|
||||
plot(lambda_nm, D_exact(1,:), 'LineWidth',2, 'DisplayName',sprintf('S_0 = 0.07'),'Color',[cols(2,:)]);
|
||||
|
||||
% plot(lambda_nm, D_lin(2,:), '-', 'LineWidth',2, 'HandleVisibility','off','Color',[cols(3,:)]);
|
||||
plot(lambda_nm, D_exact(2,:), 'LineWidth',2, 'HandleVisibility','on', 'DisplayName',sprintf('S_0 = 0.09'),'Color',[cols(4,:)]);
|
||||
plot(lambda_nm, D_exact_1(1,:), 'LineWidth',2, 'DisplayName',sprintf('$S_0=\SI{0.09}{\Dslope}$'),'Color',[0,0,0]);
|
||||
plot(lambda_nm, D_exact_1(2,:), 'LineWidth',2, 'DisplayName',sprintf('$S_0=\SI{0.09}{\Dslope}$'),'Color',[0,0,0]);
|
||||
|
||||
% plot(lambda_nm, D_exact_2(1,:), 'LineWidth',2, 'HandleVisibility','on', 'DisplayName',sprintf('$S_0=\SI{0.09}{\Dslope}$'),'Color',[0,0,0]);
|
||||
% plot(lambda_nm, D_exact_2(2,:), 'LineWidth',2, 'HandleVisibility','on', 'DisplayName',sprintf('$S_0=\SI{0.09}{\Dslope}$'),'Color',[0,0,0]);
|
||||
|
||||
xlabel('Wavelength $\lambda$ [nm]','Interpreter','latex');
|
||||
ylabel('D($\lambda$) [ps/(nm km)]','Interpreter','latex');
|
||||
@@ -44,3 +48,167 @@ xlim([min(lambda_nm) max(lambda_nm)])
|
||||
ylim([-5 5]);
|
||||
|
||||
% mat2tikz_improved('C:\Users\Silas\Documents\6971e0b65b380ca6d71c837f\02_IMDD_System\tikz\dispersion\dispersion')
|
||||
|
||||
|
||||
|
||||
%%
|
||||
|
||||
%% plot_dispersion_dual_ZDW
|
||||
% Visualizes two ZDW fiber types, each with an S0 tolerance range.
|
||||
|
||||
%% plot_dispersion_minimal_for_tikz
|
||||
% Minimal setup for clean TikZ tuning
|
||||
|
||||
% 1. Parameters
|
||||
lambda_nm = linspace(1240, 1360, 400);
|
||||
S0_range = [0.075, 0.095];
|
||||
ZDW_range = [1300, 1325];
|
||||
|
||||
% 2. Calculate Envelope and Nominal
|
||||
[S_mesh, Z_mesh] = meshgrid(S0_range, ZDW_range);
|
||||
D_all = zeros(length(lambda_nm), 4);
|
||||
for i = 1:4
|
||||
D_all(:,i) = (S_mesh(i)/4) .* (lambda_nm - (Z_mesh(i)^4)./(lambda_nm.^3));
|
||||
end
|
||||
D_env_min = min(D_all, [], 2);
|
||||
D_env_max = max(D_all, [], 2);
|
||||
D_nominal = (mean(S0_range)/4) .* (lambda_nm - (mean(ZDW_range)^4)./(lambda_nm.^3));
|
||||
|
||||
% 3. Plotting
|
||||
figure('Color','w'); hold on;
|
||||
env_col = [0.1216, 0.4706, 0.7059];
|
||||
|
||||
% Shaded area
|
||||
[hl, hp] = boundedline(lambda_nm, (D_env_min+D_env_max)/2, ...
|
||||
(D_env_max-D_env_min)/2, 'alpha', 'cmap', env_col);
|
||||
set(hl, 'YData', D_nominal, 'Color', 'k', 'LineWidth', 1.5);
|
||||
|
||||
% Generate the outline and capture the handle
|
||||
ho = outlinebounds(hl, hp);
|
||||
|
||||
% Change properties
|
||||
set(ho, 'Color', 'k', ... % Make it black
|
||||
'LineStyle', '--', ... % Make it dashed
|
||||
'LineWidth', 0.5, ... % Make it thin
|
||||
'HandleVisibility', 'off'); % Hide from legend
|
||||
|
||||
% 4. "Minimal" Measurement Lines (Tweak these in TikZ later)
|
||||
% Vertical measurement at 1290nm
|
||||
lambda_v = 1290;
|
||||
[~, idx] = min(abs(lambda_nm - lambda_v));
|
||||
y_bot = D_env_min(idx);
|
||||
y_top = D_env_max(idx);
|
||||
|
||||
% Simple line - in TikZ this will be a single \draw command
|
||||
line([lambda_v, lambda_v], [y_bot, y_top], 'Color', 'r', 'LineWidth', 1.2, 'Tag', 'VertArrow');
|
||||
text(lambda_v + 2, (y_top+y_bot)/2, '$\Delta D$', 'Color', 'r', 'Interpreter', 'latex');
|
||||
|
||||
% Horizontal measurement at D=0
|
||||
z1 = ZDW_range(1);
|
||||
z2 = ZDW_range(2);
|
||||
line([z1, z2], [0, 0], 'Color', [0.1 0.5 0.1], 'LineWidth', 1.2, 'Tag', 'HorizArrow');
|
||||
text(mean(ZDW_range), 0.5, '$\Delta \lambda_0$', 'Color', [0.1 0.5 0.1], ...
|
||||
'Interpreter', 'latex', 'HorizontalAlignment', 'center');
|
||||
|
||||
% 5. Aesthetics
|
||||
xlabel('Wavelength $\lambda$ [nm]', 'Interpreter', 'latex');
|
||||
ylabel('$D(\lambda)$ [ps/(nm km)]', 'Interpreter', 'latex');
|
||||
grid on; box on;
|
||||
xlim([1240 1360]); ylim([-5 5]);
|
||||
line(xlim, [0 0], 'Color', [0.4 0.4 0.4], 'LineStyle', '--', 'HandleVisibility', 'off');
|
||||
|
||||
% Legend using the handles we have
|
||||
legend([hp, hl], {'Worst-case Envelope', 'Nominal Realization'}, ...
|
||||
'Location', 'northwest', 'Interpreter', 'latex');
|
||||
|
||||
% To export, run:
|
||||
% matlab2tikz('dispersion.tex', 'standalone', true);
|
||||
|
||||
%% plot_dispersion_for_tikz
|
||||
% Optimized for matlab2tikz compatibility with manual arrows
|
||||
|
||||
%% plot_dispersion_statistical_envelopes
|
||||
% Visualizes hard spec limits vs. 98% statistical distribution
|
||||
|
||||
%% plot_dispersion_final_for_tikz
|
||||
lambda_nm = linspace(1240, 1360, 400);
|
||||
|
||||
% 1. Statistical & Specification Parameters
|
||||
p01_L = norminv(0.01, 1317, 2);
|
||||
p99_L = norminv(0.99, 1317, 2);
|
||||
p01_S = norminv(0.01, 0.0872, 0.0012);
|
||||
p99_S = norminv(0.99, 0.0872, 0.0012);
|
||||
|
||||
% Scenarios: [ZDW_min, ZDW_max], [S0_min, S0_max], [Color RGB]
|
||||
scenarios = { ...
|
||||
[1303, 1325], [0.075, 0.0925], [0.6510, 0.8078, 0.8902]; ... % 1. Wide Spec (Gray)
|
||||
[p01_L, p99_L], [p01_S, p99_S], [0.6, 0.6, 0.6] ... % 2. 98% Stats (Blue)
|
||||
};
|
||||
|
||||
figure('Color','w'); hold on;
|
||||
hp_handles = [];
|
||||
|
||||
% 2. Calculate and Plot Envelopes
|
||||
for k = 1:size(scenarios, 1)
|
||||
Zr = scenarios{k,1};
|
||||
Sr = scenarios{k,2};
|
||||
col = scenarios{k,3};
|
||||
|
||||
[S_mesh, Z_mesh] = meshgrid(Sr, Zr);
|
||||
D_all = zeros(length(lambda_nm), 4);
|
||||
for i = 1:4
|
||||
D_all(:,i) = (S_mesh(i)/4) .* (lambda_nm - (Z_mesh(i)^4)./(lambda_nm.^3));
|
||||
end
|
||||
|
||||
D_min_env = min(D_all, [], 2);
|
||||
D_max_env = max(D_all, [], 2);
|
||||
|
||||
[hl, hp] = boundedline(lambda_nm, (D_min_env+D_max_env)/2, (D_max_env-D_min_env)/2, ...
|
||||
'cmap','alpha', col);
|
||||
|
||||
hp_handles(k) = hp;
|
||||
set(hl, 'Visible', 'off');
|
||||
|
||||
|
||||
if k == 1
|
||||
D_wide_min = D_min_env;
|
||||
D_wide_max = D_max_env;
|
||||
% ho = outlinebounds(hl, hp);
|
||||
% % Change properties
|
||||
% set(ho, 'Color', 'k', ... % Make it black
|
||||
% 'LineStyle', '--', ... % Make it dashed
|
||||
% 'LineWidth', 1, ... % Make it thin
|
||||
% 'HandleVisibility', 'off'); % Hide from legend
|
||||
% Capture Wide Spec (k=1) bounds for the TikZ measurement lines
|
||||
hp.FaceAlpha = 0.5;
|
||||
else
|
||||
hp.FaceAlpha = 0.8;
|
||||
end
|
||||
end
|
||||
|
||||
% 3. Nominal Line
|
||||
D_nom = (0.0872/4) .* (lambda_nm - (1317^4)./(lambda_nm.^3));
|
||||
h_nom = plot(lambda_nm, D_nom, 'k', 'LineWidth', 1,'LineStyle','-');
|
||||
|
||||
% 4. Minimal Lines for TikZ (Measuring Wide Spec)
|
||||
lambda_v = 1290;
|
||||
[~, idx_v] = min(abs(lambda_nm - lambda_v));
|
||||
% Vertical line showing full Wide Spec dispersion range at 1290nm
|
||||
line([lambda_v, lambda_v], [D_wide_min(idx_v), D_wide_max(idx_v)], 'Color', 'k', 'Tag', 'VertArrow','LineWidth', 1);
|
||||
% Horizontal line showing Wide Spec ZDW range at D=0
|
||||
% line([1303, 1325], [0, 0], 'Color', 'k', 'Tag', 'HorizArrow','LineWidth', 1);
|
||||
|
||||
% 5. Aesthetics & Legend
|
||||
xlabel('Wavelength $\lambda$ [nm]', 'Interpreter', 'latex');
|
||||
ylabel('$D(\lambda)$ [ps/(nm km)]', 'Interpreter', 'latex');
|
||||
grid on; box on;
|
||||
xlim([1240 1360]); ylim([-5 5]);
|
||||
|
||||
leg_labels = { ...
|
||||
'$\lambda_0 \in [1303, 1325], S_0 \in [0.075, 0.0925]$', ...
|
||||
'$\lambda_0 \in [1312, 1322], S_0 \in [0.084, 0.090]$', ...
|
||||
'$\lambda_0 = 1317, S_0 = 0.0872$'};
|
||||
legend([hp_handles, h_nom], leg_labels, 'Location', 'northwest', 'Interpreter', 'latex', 'FontSize', 8);
|
||||
|
||||
% Export command (uncomment to use)
|
||||
% mat2tikz_improved('C:\Users\Silas\Documents\6971e0b65b380ca6d71c837f\02_IMDD_System\tikz\dispersion\dispersion_slope.tikz')
|
||||
@@ -0,0 +1,183 @@
|
||||
%% ------------------------------------------------------------
|
||||
% Contour plot: λ_null as function of bandwidth (f_target) and reach (L)
|
||||
% ------------------------------------------------------------
|
||||
|
||||
% Parameters
|
||||
lambda0 = 1310e-9; % [m]
|
||||
S0 = 0.09; % [ps/(nm²·km)]
|
||||
c = physconst('lightspeed');
|
||||
|
||||
% Sweep dimensions
|
||||
f_targets = linspace(50e9, 130e9, 200); % [Hz] (x-axis)
|
||||
L_values = linspace(0.5e3, 15e3, 200); % [m] (y-axis)
|
||||
|
||||
lambda_surface = zeros(numel(L_values), numel(f_targets));
|
||||
Dacc_surface = zeros(numel(L_values), numel(f_targets));
|
||||
|
||||
% Outer loop over fiber length (since L must be scalar)
|
||||
for iL = 1:numel(L_values)
|
||||
L = L_values(iL);
|
||||
[lambda_vec, Dacc_vec] = lambda_for_first_null_full(f_targets, L, lambda0, S0);
|
||||
|
||||
% Store the 1x absolute offset |lambda - lambda0|
|
||||
lambda_surface(iL, :) = abs(lambda0 - lambda_vec);
|
||||
Dacc_surface(iL, :) = Dacc_vec;
|
||||
end
|
||||
|
||||
% Convert for plotting
|
||||
lambda_surface_nm = lambda_surface * 1e9; % [nm]
|
||||
L_km = L_values / 1000; % [km]
|
||||
f_GHz = f_targets / 1e9; % [GHz]
|
||||
|
||||
%% Contour plot
|
||||
figure('Color','w');
|
||||
hold on;
|
||||
|
||||
% Define wavelength contour levels [nm]
|
||||
% Focus on a clean range of 1x offset values
|
||||
lambda_levels = unique([50:-10:30, 30:-5:5]);
|
||||
|
||||
% Contour plot
|
||||
[C,h] = contour(f_GHz, L_km, lambda_surface_nm, lambda_levels, ...
|
||||
'LineWidth', 1.2, ...
|
||||
'ShowText', 'off');
|
||||
|
||||
% Colormap: Modern Blue palette with light colors removed for visibility
|
||||
cmap_full = cbrewer2('Blues', 40);
|
||||
colormap(cmap_full(10:end, :));
|
||||
clim([min(lambda_levels) max(lambda_levels)]);
|
||||
cb = colorbar;
|
||||
ylabel(cb, '$\Delta \lambda$ [nm]', 'Interpreter', 'latex');
|
||||
|
||||
% --- MANUAL TEXTBOX ANNOTATIONS ---
|
||||
% Find placement along the first-null curve for each level
|
||||
for i = 1:length(lambda_levels)
|
||||
lvl = lambda_levels(i);
|
||||
|
||||
% Re-calculate the specific (f, L) curve for this delta-lambda
|
||||
lambda_target = lambda0 - (lvl * 1e-9);
|
||||
LHS = -( (S0*1e3) / 4 ) * (lambda_target - (lambda0^4)/(lambda_target^3)) * lambda_target^2;
|
||||
const_val = (c*0.5) / LHS;
|
||||
|
||||
f_curve_GHz = linspace(min(f_GHz), max(f_GHz), 500);
|
||||
L_curve_km = const_val ./ (f_curve_GHz * 1e9).^2 / 1000;
|
||||
|
||||
% Filter for points within the plot axes
|
||||
in_bounds = find(L_curve_km >= min(L_km)*1.1 & L_curve_km <= max(L_km)*0.9 & ...
|
||||
f_curve_GHz >= min(f_GHz)*1.1 & f_curve_GHz <= max(f_GHz)*0.9);
|
||||
|
||||
if ~isempty(in_bounds)
|
||||
% Specific alternating pattern for weight to minimize overlapping
|
||||
if i < 9
|
||||
weight = 0.05;
|
||||
else
|
||||
weight = 0.05 + 0.1 * mod(i, 2);
|
||||
end
|
||||
idx = in_bounds(max(1, min(length(in_bounds), round(length(in_bounds) * weight))));
|
||||
|
||||
text(f_curve_GHz(idx), L_curve_km(idx), sprintf('%g nm', lvl), ...
|
||||
'Color', 'k', 'BackgroundColor', 'w', 'Margin', 1.5, ...
|
||||
'HorizontalAlignment', 'center', 'VerticalAlignment', 'middle', ...
|
||||
'EdgeColor', 'k', 'FontSize', 9);
|
||||
end
|
||||
end
|
||||
|
||||
% Axis formatting
|
||||
xlabel('Signal Bandwidth [GHz]', 'FontSize', 11);
|
||||
ylabel('Fiber length [km]', 'FontSize', 11);
|
||||
xticks(min(f_GHz):10:max(f_GHz));
|
||||
yticks(min(L_km):2.5:max(L_km));
|
||||
grid on; box on;
|
||||
axis([min(f_GHz) max(f_GHz) min(L_km) max(L_km)]);
|
||||
|
||||
|
||||
%% Optional: overlay accumulated-dispersion contours
|
||||
if 0
|
||||
hold on;
|
||||
min_D = min(Dacc_surface(:), [], 'omitnan');
|
||||
max_D = max(Dacc_surface(:), [], 'omitnan');
|
||||
% Calculate 3 integer levels well within the data range
|
||||
D_levels = unique(round(linspace(min_D*0.8, max_D*0.8, 3)));
|
||||
|
||||
[CS, h] = contour(f_GHz, L_km, Dacc_surface, D_levels, 'k--', 'LineWidth', 0.8);
|
||||
clabel(CS, h, 'Color','k', 'FontSize',8);
|
||||
end
|
||||
|
||||
|
||||
%% Export
|
||||
% Hier erzwingen wir die rote Colormap für pgfplots, damit mat2tikz es nicht blau exportiert!
|
||||
mat2tikz_improved("C:/Users/Silas/Documents/6971e0b65b380ca6d71c837f/02_IMDD_System/tikz/dispersion/dispersion_power_fading_contour2.tikz");
|
||||
|
||||
function [lambda_vec, Dacc_vec] = lambda_for_first_null_full(f_target, L, lambda0, S0)
|
||||
% lambda_for_first_null_full (stable, single-branch + validity checks)
|
||||
% --------------------------------------------------------------------
|
||||
% Computes the wavelength(s) at which the first IM/DD fading null
|
||||
% occurs at frequency/ies f_target using the full dispersion model:
|
||||
%
|
||||
% D(lambda) = (S0/4)*(lambda - lambda0^4 / lambda^3)
|
||||
%
|
||||
% Restricted to the NORMAL-dispersion branch (λ < λ0),
|
||||
% and valid only in the O-band (1260–1360 nm).
|
||||
%
|
||||
% Inputs:
|
||||
% f_target - scalar or vector of target null frequencies [Hz]
|
||||
% L - fiber length [m]
|
||||
% lambda0 - zero-dispersion wavelength (ZDW) [m]
|
||||
% S0 - dispersion slope at ZDW [ps/(nm²·km)]
|
||||
%
|
||||
% Outputs:
|
||||
% lambda_vec - wavelength(s) [m] where first null occurs (clamped to O-band)
|
||||
% Dacc_vec - accumulated dispersion(s) [ps/nm] (NaN if out of valid range)
|
||||
% --------------------------------------------------------------------
|
||||
|
||||
c = physconst('lightspeed');
|
||||
S0_si = S0 * 1e3; % ps/(nm²·km) -> s/(m³)
|
||||
|
||||
% Define O-band boundaries (in meters)
|
||||
lambda_min = 1255e-9;
|
||||
lambda_max = 1361e-9;
|
||||
|
||||
% Force column vector
|
||||
f_target = f_target(:);
|
||||
N = numel(f_target);
|
||||
|
||||
lambda_vec = NaN(N,1);
|
||||
Dacc_vec = NaN(N,1);
|
||||
|
||||
|
||||
|
||||
for k = 1:N
|
||||
RHS = c * 0.5 / (f_target(k)^2 * L);
|
||||
|
||||
% Normal-dispersion branch (λ < λ0)
|
||||
fun = @(lambda) -(S0_si/4).*(lambda - (lambda0^4)./(lambda.^3)).*lambda.^2 - RHS;
|
||||
|
||||
% Limit the search to [λ_min, λ0)
|
||||
try
|
||||
lambda_sol = fzero(fun, [lambda_min, lambda0 * 0.999]);
|
||||
catch
|
||||
% If the zero is not within bounds, skip this point
|
||||
lambda_sol = NaN;
|
||||
end
|
||||
|
||||
% Validate solution
|
||||
if isnan(lambda_sol) || lambda_sol < lambda_min || lambda_sol > lambda_max
|
||||
lambda_vec(k) = NaN;
|
||||
Dacc_vec(k) = NaN;
|
||||
continue
|
||||
end
|
||||
|
||||
% Compute D(lambda) and accumulated dispersion
|
||||
D_lambda = (S0_si/4) * (lambda_sol - (lambda0^4)/(lambda_sol^3)) / 1e-6; % ps/(nm·km)
|
||||
Dacc_val = D_lambda * (L/1000); % ps/nm
|
||||
|
||||
% Sanity bound on dispersion (avoid unphysical > ±100 ps/nm)
|
||||
if abs(Dacc_val) > 100
|
||||
lambda_vec(k) = NaN;
|
||||
Dacc_vec(k) = NaN;
|
||||
else
|
||||
lambda_vec(k) = lambda_sol;
|
||||
Dacc_vec(k) = Dacc_val;
|
||||
end
|
||||
end
|
||||
end
|
||||
142
Functions/Theory/Dissertation/dispersion_contour_bivariate.m
Normal file
142
Functions/Theory/Dissertation/dispersion_contour_bivariate.m
Normal file
@@ -0,0 +1,142 @@
|
||||
% Festen Betriebsparameter
|
||||
lambda = 1290; % nm
|
||||
L = 1; % km
|
||||
mu_zwd = 1317; % nm
|
||||
sigma_zwd = 2; % nm
|
||||
mu_s0 = 0.0872; % ps / nm2 km
|
||||
sigma_s0 = 0.0012; % ps / nm2 km
|
||||
rho = -0.5; % Korrelation
|
||||
|
||||
% Gitter für lambda0 und S0
|
||||
lambda0_vec = linspace(mu_zwd-10, mu_zwd+10, 100);
|
||||
S0_vec = linspace(mu_s0-0.01, mu_s0+0.01, 100);
|
||||
|
||||
% Korrigierte meshgrid Reihenfolge
|
||||
[S0, Lambda0] = meshgrid(S0_vec, lambda0_vec);
|
||||
|
||||
% Dispersion berechnen
|
||||
D = (S0./4) .* ( lambda - (Lambda0.^4)./(lambda^3) ) * L;
|
||||
|
||||
% 2D Verteilung (Bivariate Gauss) berechnen
|
||||
Z_x = (S0 - mu_s0) / sigma_s0;
|
||||
Z_y = (Lambda0 - mu_zwd) / sigma_zwd;
|
||||
PDF_2D = exp(-1 / (2 * (1 - rho^2)) * (Z_x.^2 - 2 * rho .* Z_x .* Z_y + Z_y.^2));
|
||||
|
||||
%% Plot zusammenbauen
|
||||
figure('Color','w');
|
||||
hold on % EINZIGES hold on für den gesamten Plot!
|
||||
|
||||
% --- 1. ZUERST: 2D Verteilung (Bivariate Gauss) "ganz unten" ---
|
||||
numLevels_2D = 6;
|
||||
colors_2D = cbrewer2('Greys', numLevels_2D+0);
|
||||
colors_2D = colors_2D(1:numLevels_2D,:);
|
||||
% Nur die Anzahl der Level übergeben!
|
||||
[C2, h2] = contourf(S0, Lambda0, PDF_2D, numLevels_2D);
|
||||
h2.HandleVisibility='off';
|
||||
% Colormap für die rote Fläche setzen
|
||||
|
||||
colormap(gcf, colors_2D(1:end,:));
|
||||
try
|
||||
clim([min(PDF_2D(:)), max(PDF_2D(:))]);
|
||||
catch
|
||||
caxis([min(PDF_2D(:)), max(PDF_2D(:))]);
|
||||
end
|
||||
|
||||
h2.EdgeColor = 'none'; % Keine schwarzen Ränder
|
||||
|
||||
% --- 2. DARÜBER: Dispersions-Konturlinien ---
|
||||
numLevels_D = 9;
|
||||
% Erzeuge glatte, auf 1 Nachkommastelle gerundete Werte
|
||||
levels_D = round(linspace(min(D(:)), max(D(:)), numLevels_D), 1);
|
||||
levels_D = unique(levels_D);
|
||||
cmap_bg = flip(cbrewer2('Blues', length(levels_D)+3));
|
||||
|
||||
for i = 1:length(levels_D)
|
||||
% Konturlinien zeichnen (explizit Schwarz)
|
||||
[C,h] = contour(S0, Lambda0, D, [levels_D(i), levels_D(i)], ...
|
||||
'Color', cmap_bg(i,:),...
|
||||
'LineWidth', 1, ...
|
||||
'ShowText', 'off','handlevisibility','off');
|
||||
if i == 1
|
||||
h.HandleVisibility='on';
|
||||
h.DisplayName='$D(\lambda=1290)$';
|
||||
end
|
||||
end
|
||||
|
||||
% --- 3. MANUELLE TEXTBOXEN AUF DEN LINIEN ---
|
||||
% Wähle eine feste S0-Position für alle Beschriftungen (z.B. bei 0.082)
|
||||
S0_label = 0.095;
|
||||
|
||||
for i = 1:length(levels_D)
|
||||
% Berechne exakte ZDW (Y-Koordinate) durch Umstellen der D-Formel:
|
||||
% Lambda0 = (lambda^3 * (lambda - 4*D / (S0 * L)))^(1/4)
|
||||
zdw_label = (lambda^3 * (lambda - 4*levels_D(i) / (S0_label * L)))^0.25;
|
||||
|
||||
% Nur zeichnen, wenn der Punkt auch im sichtbaren Plot-Bereich liegt
|
||||
if zdw_label >= min(lambda0_vec) && zdw_label <= max(lambda0_vec)
|
||||
text(S0_label, zdw_label, sprintf('%0.1f', levels_D(i)), ...
|
||||
'Color', 'k', ...
|
||||
'BackgroundColor', 'w', ... % Weiße Box überdeckt die schwarze Linie!
|
||||
'Margin', 2, ... % Abstand der Box zum Text
|
||||
'HorizontalAlignment', 'center', ...
|
||||
'VerticalAlignment', 'middle', ...
|
||||
'FontSize', 10);
|
||||
end
|
||||
end
|
||||
|
||||
% --- 3. GANZ OBEN: Randverteilungen (1D Gauss) an den Achsen ---
|
||||
% S0 Verteilung (unten)
|
||||
s0_vals = linspace(min(S0_vec), max(S0_vec), 500);
|
||||
gauss_s0 = exp(-0.5*((s0_vals - mu_s0)/sigma_s0).^2);
|
||||
scale_s0 = 4;
|
||||
y_s0_base = min(lambda0_vec);
|
||||
plot(s0_vals, y_s0_base + gauss_s0 * scale_s0, 'LineWidth', 1,'LineStyle','--','DisplayName','$S_0$','Color',[0,0,0],'HandleVisibility','off');
|
||||
% ANNOTATION S0: Automatisch platziert leicht über dem Peak
|
||||
str_s0 = sprintf('\\mu_{S0} = %.4f\n\\sigma_{S0} = %.4f', mu_s0, sigma_s0);
|
||||
text(0.0915,1308, str_s0, ...
|
||||
'Interpreter', 'tex', ...
|
||||
'HorizontalAlignment', 'left', ... % Entspricht 'right' in TikZ
|
||||
'VerticalAlignment', 'middle', ...
|
||||
'BackgroundColor', 'w', ... % Entspricht 'fill=white'
|
||||
'EdgeColor', 'k', ... % Entspricht 'draw=black'
|
||||
'Margin', 1, ... % Entspricht 'inner sep=1pt'
|
||||
'FontSize', 10);
|
||||
|
||||
% ZDW Verteilung (links)
|
||||
% ZDW Verteilung (links)
|
||||
zwd_vals = linspace(min(lambda0_vec), max(lambda0_vec), 500);
|
||||
gauss_zwd = exp(-0.5*((zwd_vals - mu_zwd)/sigma_zwd).^2);
|
||||
scale_zwd = 0.005;
|
||||
x_zwd_base = min(S0_vec);
|
||||
plot(x_zwd_base + gauss_zwd * scale_zwd, zwd_vals, 'LineWidth', 1,'LineStyle','--','DisplayName','ZDW','Color',[0,0,0],'HandleVisibility','off');
|
||||
|
||||
% ANNOTATION ZDW: Automatisch platziert leicht rechts neben dem Peak
|
||||
str_zwd = sprintf('\\mu_{ZDW} = %.1f\n\\sigma_{ZDW} = %.1f', mu_zwd, sigma_zwd);
|
||||
text(0.079,1312, str_zwd, ...
|
||||
'Interpreter', 'tex', ...
|
||||
'HorizontalAlignment', 'left', ... % Entspricht 'right' in TikZ
|
||||
'VerticalAlignment', 'middle', ...
|
||||
'BackgroundColor', 'w', ... % Entspricht 'fill=white'
|
||||
'EdgeColor', 'k', ... % Entspricht 'draw=black'
|
||||
'Margin', 1, ... % Entspricht 'inner sep=1pt'
|
||||
'FontSize', 10);
|
||||
|
||||
% Hilfslinien für die Mittelwerte
|
||||
% xline(mu_s0,'LineWidth',0.5,'HandleVisibility','off','LineStyle',':');
|
||||
% yline(mu_zwd,'LineWidth',0.5,'HandleVisibility','off','LineStyle',':');
|
||||
|
||||
% --- Achsenbeschriftung, Titel & Formatierung ---
|
||||
xlabel('$S_0$ [$\nicefrac{\text{ps}}{(\text{nm}^2\text{ km})}$]', 'FontSize', 12);
|
||||
ylabel('ZDW [nm]', 'FontSize', 12);
|
||||
% title(sprintf('Dispersion: %d km; %d nm', L, lambda), 'FontSize', 14);
|
||||
|
||||
grid on
|
||||
% Exakte Begrenzung, damit die Randverteilungen bündig anliegen
|
||||
axis([min(S0_vec) max(S0_vec) min(lambda0_vec) max(lambda0_vec)]);
|
||||
|
||||
% legend;
|
||||
hold off % EINZIGES hold off ganz am Ende!
|
||||
|
||||
%% Export
|
||||
% Hier erzwingen wir die rote Colormap für pgfplots, damit mat2tikz es nicht blau exportiert!
|
||||
mat2tikz_improved("C:/Users/Silas/Documents/6971e0b65b380ca6d71c837f/02_IMDD_System/tikz/dispersion/dispersion_contour_bi2.tikz");
|
||||
@@ -1,124 +0,0 @@
|
||||
%% ============================================================
|
||||
% SNR vs Optical Input Power (dBm) — Shot vs Thermal vs Combined
|
||||
% - Generate optical field as sine with target RMS power (verified)
|
||||
% - Magnitude-square detection -> optical power
|
||||
% - Photocurrent = Rd * P
|
||||
% - Add shot noise + thermal noise (white, PSD-based)
|
||||
% - Compute SNR for: shot-only, thermal-only, combined
|
||||
% ============================================================
|
||||
|
||||
|
||||
% Constants
|
||||
k = Constant.Boltzmann;
|
||||
q = Constant.ElementaryCharge;
|
||||
|
||||
% Receiver / PD parameters
|
||||
T = 20 + 273.15; % K
|
||||
R = 50; % Ohm (front-end/load)
|
||||
Rd = 0.7; % A/W (responsivity)
|
||||
Be = 100e9; % Hz (electrical noise bandwidth)
|
||||
Id = 0; % A (dark current, optional)
|
||||
|
||||
% Sampling for time-domain demo (needs to be >> Be)
|
||||
fs = 1e12; % Hz
|
||||
N = 2^14; % samples
|
||||
t = (0:N-1).'/fs;
|
||||
|
||||
% Choose an electrical tone within bandwidth (arbitrary for demo)
|
||||
f0 = 10e9; % Hz
|
||||
|
||||
% Thermal current PSD (two-sided) and variance in Be
|
||||
Si_th = 4*k*T/R; % A^2/Hz (two-sided)
|
||||
sigma2_th = Si_th * Be; % A^2
|
||||
sigma_th = sqrt(sigma2_th); % A_rms
|
||||
|
||||
% Optical input power sweep (in dBm)
|
||||
P_dBm = linspace(-40, 10, 300);
|
||||
P_W = 10.^((P_dBm - 30)/10); % W
|
||||
|
||||
% Pre-allocate results
|
||||
SNR_shot_dB = zeros(size(P_dBm));
|
||||
SNR_th_dB = zeros(size(P_dBm));
|
||||
SNR_tot_dB = zeros(size(P_dBm));
|
||||
|
||||
% --- Main loop over optical input power
|
||||
for ii = 1:numel(P_W)
|
||||
Pavg = P_W(ii);
|
||||
|
||||
% Optical field with RMS power = Pavg:
|
||||
% Let x(t) be the optical field amplitude such that |x|^2 has mean Pavg.
|
||||
% Use a sinusoid: x(t) = A*sin(2*pi*f0*t), then mean(|x|^2)=A^2/2.
|
||||
A = sqrt(2*Pavg); % -> mean(|x|^2) = Pavg
|
||||
|
||||
x = A * sin(2*pi*f0*t); % "optical field" (real for simplicity)
|
||||
|
||||
% Verify RMS/mean power numerically (optional)
|
||||
P_meas = mean(abs(x).^2); % should be ~ Pavg
|
||||
|
||||
a = P_meas - Pavg;
|
||||
assert(a<1);
|
||||
|
||||
% Magnitude-square detection -> optical power waveform
|
||||
Popt = abs(x).^2; % W (instantaneous)
|
||||
|
||||
% Photocurrent waveform (includes DC + 2f0 component for this demo)
|
||||
I_sig = Rd * Popt; % A
|
||||
|
||||
% Define "signal power" as the mean-squared photocurrent due to signal
|
||||
% (for this sine-squared waveform, it's OK for a demo SNR definition)
|
||||
P_sig = mean(I_sig.^2);
|
||||
|
||||
% -------- Shot noise (white) --------
|
||||
% Two-sided PSD: Si_shot = 2*q*(I_photo + I_dark)
|
||||
% For this demo, use average current to set the white-noise level:
|
||||
Ibar = mean(I_sig) + Id;
|
||||
Si_shot = 2*q*Ibar; % A^2/Hz
|
||||
sigma2_sh = Si_shot * Be; % A^2
|
||||
sigma_sh = sqrt(sigma2_sh); % A_rms
|
||||
|
||||
n_sh = sigma_sh * randn(N,1); % time-domain shot noise
|
||||
|
||||
% -------- Thermal noise (white Gaussian) --------
|
||||
n_th = sigma_th * randn(N,1); % time-domain thermal noise
|
||||
|
||||
% Noise powers (mean-square) in time domain
|
||||
Pn_sh = mean(n_sh.^2);
|
||||
Pn_th = mean(n_th.^2);
|
||||
Pn_tot = mean((n_sh + n_th).^2); % ~ Pn_sh + Pn_th (independent)
|
||||
|
||||
% SNRs
|
||||
SNR_shot = P_sig / Pn_sh;
|
||||
SNR_th = P_sig / Pn_th;
|
||||
SNR_tot = P_sig / Pn_tot;
|
||||
|
||||
SNR_shot_dB(ii) = 10*log10(SNR_shot);
|
||||
SNR_th_dB(ii) = 10*log10(SNR_th);
|
||||
SNR_tot_dB(ii) = 10*log10(SNR_tot);
|
||||
|
||||
% Optional sanity check (can comment out)
|
||||
% if ii == 1
|
||||
% fprintf('Check Pavg target/meas: %.3g W / %.3g W\n', Pavg, P_meas);
|
||||
% end
|
||||
end
|
||||
|
||||
% Plot comparison
|
||||
figure(1); clf; hold on;
|
||||
plot(P_dBm, SNR_shot_dB, 'LineWidth', 1.5);
|
||||
plot(P_dBm, SNR_th_dB, 'LineWidth', 1.5);
|
||||
plot(P_dBm, SNR_tot_dB, 'LineWidth', 1.5);
|
||||
grid on;
|
||||
xlabel('Optical input power [dBm]');
|
||||
ylabel('SNR [dB]');
|
||||
title('SNR vs Optical Input Power: Shot vs Thermal vs Combined');
|
||||
legend('Shot-noise only','Thermal-noise only','Shot + Thermal','Location','best');
|
||||
|
||||
% Print example at 0 dBm
|
||||
[~,idx0] = min(abs(P_dBm - 0));
|
||||
fprintf('At 0 dBm:\n');
|
||||
fprintf(' SNR (shot only) = %.2f dB\n', SNR_shot_dB(idx0));
|
||||
fprintf(' SNR (thermal only) = %.2f dB\n', SNR_th_dB(idx0));
|
||||
fprintf(' SNR (combined) = %.2f dB\n', SNR_tot_dB(idx0));
|
||||
|
||||
% Also print NEP based on thermal PSD (constant NEP_th)
|
||||
NEP_th = sqrt(Si_th)/Rd; % W/sqrt(Hz)
|
||||
fprintf('Thermal NEP = %.3g W/sqrt(Hz)\n', NEP_th);
|
||||
@@ -17,8 +17,8 @@ D_si = D_lambda * 1e-6; % → s/m²
|
||||
b2 = -D_si * lambda^2 / (2*pi*c); % s²/m
|
||||
|
||||
%% Frequency grid
|
||||
f_max = 150e9;
|
||||
f = linspace(0, f_max, 4000); % [Hz]
|
||||
f_max = 200e9;
|
||||
f = linspace(0, f_max, 5000); % [Hz]
|
||||
|
||||
%% IM/DD transfer function (power fading)
|
||||
phi = 2*pi^2 * b2 * f.^2 * L;
|
||||
@@ -26,17 +26,22 @@ H = abs(cos(phi));
|
||||
|
||||
%% Plot
|
||||
figure('Color','w');
|
||||
plot(f/1e9, 10*log10(H), 'LineWidth', 1.8,'Color','black');
|
||||
plot(f/1e9, 10*log10(H), 'LineWidth', 1,'Color','black');
|
||||
grid on; box on;
|
||||
xlabel('Frequency [GHz]');
|
||||
ylabel('Magnitude [dB]');
|
||||
% title(sprintf('IM/DD Power Fading: 10 km; 1275nm', lambda*1e9, L/1000),"Interpreter","latex");
|
||||
ylim([-30 0]);
|
||||
ylim([-20 0]);
|
||||
|
||||
%% Mark analytic first-null frequency
|
||||
f_null = sqrt(c*(0.5)/(abs(D_si)*lambda^2*L));
|
||||
xline(f_null/1e9, 'r--', 'LineWidth', 1.2, ...
|
||||
'Label', sprintf('f_{null}=%.1f GHz', f_null/1e9), ...
|
||||
'LabelOrientation', 'horizontal', 'LabelVerticalAlignment', 'bottom');
|
||||
%% Mark analytic null frequencies up to order 5
|
||||
max_order = 3;
|
||||
for n = 0:max_order
|
||||
f_null_n = sqrt( c*(2*n + 1)/(2*abs(D_si)*lambda^2*L) );
|
||||
xline(f_null_n/1e9, '--', 'LineWidth', 1.2, ...
|
||||
'Color', [0.1216, 0.4706, 0.7059]);
|
||||
text(f_null_n/1e9, -15, sprintf('$f_{\\mathrm{null}, %d}=%.1f$ GHz', n, f_null_n/1e9), ...
|
||||
'BackgroundColor', 'w', 'EdgeColor', 'k', 'Interpreter', 'latex', ...
|
||||
'HorizontalAlignment', 'center', 'VerticalAlignment', 'middle');
|
||||
end
|
||||
|
||||
mat2tikz_improved('C:\Users\Silas\Documents\6971e0b65b380ca6d71c837f\02_IMDD_System\tikz\dispersion\power_fading.tikz')
|
||||
% mat2tikz_improved('C:\Users\Silas\Documents\6971e0b65b380ca6d71c837f\02_IMDD_System\tikz\dispersion\power_fading.tikz')
|
||||
@@ -1,68 +0,0 @@
|
||||
%% ============================================================
|
||||
% IM/DD Power Fading Evolution GIF (1 km -> 20 km)
|
||||
% Uses the provided GifWriter (serial mode)
|
||||
% ============================================================
|
||||
|
||||
clear; close all; clc;
|
||||
|
||||
%% Fiber and system parameters
|
||||
lambda0 = 1310e-9; % zero-dispersion wavelength [m]
|
||||
lambda = 1275e-9; % operating wavelength [m]
|
||||
S0 = 0.09; % dispersion slope [ps/(nm^2·km)]
|
||||
c = physconst('lightspeed');
|
||||
|
||||
%% Derived quantities (length-independent)
|
||||
D_lambda = (S0/4) * (lambda*1e9 - (lambda0*1e9)^4/(lambda*1e9)^3); % ps/(nm·km)
|
||||
D_si = D_lambda * 1e-6; % s/m^2
|
||||
b2 = -D_si * lambda^2 / (2*pi*c); % s^2/m
|
||||
|
||||
%% Frequency grid
|
||||
f_max = 150e9;
|
||||
f = linspace(0, f_max, 4000); % [Hz]
|
||||
|
||||
%% Figure setup (keep it stable for nicer GIFs)
|
||||
fig = figure('Color','w');
|
||||
ax = axes(fig); %#ok<LAXES>
|
||||
hold(ax,'on'); grid(ax,'on'); box(ax,'on');
|
||||
xlabel(ax,'Frequency [GHz]');
|
||||
ylabel(ax,'Magnitude [dB]');
|
||||
ylim(ax,[-30 0]);
|
||||
xlim(ax,[0 f_max/1e9]);
|
||||
|
||||
%% GIF writer (serial mode; simplest)
|
||||
g = GifWriter('Name','power_fading_evolution', 'DelayTime',0.12, 'Parallel',false);
|
||||
|
||||
%% Loop: 1 km to 20 km
|
||||
L = [1:20,19:-1:1];
|
||||
for L_km = L
|
||||
L_meter = L_km * 1e3; % [m]
|
||||
|
||||
% IM/DD transfer function (power fading)
|
||||
phi = 2*pi^2 * b2 * f.^2 * L_meter;
|
||||
H = abs(cos(phi));
|
||||
HdB = 10*log10(max(H, 1e-12)); % avoid -Inf for deep notches
|
||||
|
||||
% Clear and redraw (stable axes)
|
||||
cla(ax);
|
||||
|
||||
plot(ax, f/1e9, HdB, 'LineWidth', 1.8, 'Color','black');
|
||||
|
||||
% Analytic first-null frequency marker
|
||||
f_null = sqrt(c*(0.5)/(abs(D_si)*lambda^2*L_meter));
|
||||
xline(ax, f_null/1e9, 'r--', 'LineWidth', 1.2, ...
|
||||
'Label', sprintf('f_{null}=%.1f GHz', f_null/1e9), ...
|
||||
'LabelOrientation','horizontal', ...
|
||||
'LabelVerticalAlignment','bottom');
|
||||
|
||||
title(ax, sprintf('Power Fading for: %.0f km @ 1275 nm', L_km));
|
||||
|
||||
drawnow;
|
||||
|
||||
% Add frame to GIF
|
||||
g.addFrame(fig);
|
||||
end
|
||||
|
||||
%% Done
|
||||
g.compile(fig.Number);
|
||||
|
||||
disp(fullfile(g.OutputDir, sprintf('%s_fig_%d.gif', g.Name, fig.Number)));
|
||||
@@ -1,69 +0,0 @@
|
||||
%% ============================================================
|
||||
% IM/DD Power Fading Evolution vs Wavelength (L = 10 km)
|
||||
% ============================================================
|
||||
|
||||
clear; close all; clc;
|
||||
|
||||
%% Fixed fiber parameters
|
||||
lambda0 = 1310e-9; % zero-dispersion wavelength [m]
|
||||
S0 = 0.09; % dispersion slope [ps/(nm^2·km)]
|
||||
L = 10e3; % fiber length FIXED [m]
|
||||
c = physconst('lightspeed');
|
||||
|
||||
%% Frequency grid
|
||||
f_max = 150e9;
|
||||
f = linspace(0, f_max, 4000); % [Hz]
|
||||
|
||||
%% Figure setup (stable axes for clean GIF)
|
||||
fig = figure('Color','w');
|
||||
ax = axes(fig);
|
||||
hold(ax,'on'); grid(ax,'on'); box(ax,'on');
|
||||
xlabel(ax,'Frequency [GHz]');
|
||||
ylabel(ax,'Magnitude [dB]');
|
||||
ylim(ax,[-30 0]);
|
||||
xlim(ax,[0 f_max/1e9]);
|
||||
|
||||
%% GIF writer
|
||||
g = GifWriter('Name','power_fading_vs_wavelength', ...
|
||||
'DelayTime',0.12, ...
|
||||
'Parallel',false);
|
||||
|
||||
%% Wavelength sweep (around ZDW)
|
||||
lambda_vec = linspace(1260e-9, 1360e-9, 25); % 1260–1360 nm
|
||||
|
||||
for k = 1:length(lambda_vec)
|
||||
|
||||
lambda = lambda_vec(k);
|
||||
|
||||
%% Dispersion for current wavelength
|
||||
D_lambda = (S0/4) * (lambda*1e9 - (lambda0*1e9)^4/(lambda*1e9)^3); % ps/(nm·km)
|
||||
D_si = D_lambda * 1e-6; % s/m^2
|
||||
b2 = -D_si * lambda^2 / (2*pi*c); % s^2/m
|
||||
|
||||
%% Power fading transfer function
|
||||
phi = 2*pi^2 * b2 * f.^2 * L;
|
||||
H = abs(cos(phi));
|
||||
HdB = 10*log10(max(H, 1e-12));
|
||||
|
||||
cla(ax)
|
||||
plot(ax, f/1e9, HdB, 'LineWidth',1.8,'Color','black');
|
||||
|
||||
%% First-null frequency
|
||||
if abs(D_si) > 0
|
||||
f_null = sqrt(c*(0.5)/(abs(D_si)*lambda^2*L));
|
||||
xline(ax, f_null/1e9, 'r--', 'LineWidth',1.2, ...
|
||||
'Label', sprintf('f_{null}=%.1f GHz', f_null/1e9), ...
|
||||
'LabelOrientation','horizontal', ...
|
||||
'LabelVerticalAlignment','bottom');
|
||||
end
|
||||
|
||||
title(ax, sprintf('Power Fading for: 10 km @ %.0f nm', lambda*1e9));
|
||||
|
||||
drawnow;
|
||||
g.addFrame(fig);
|
||||
end
|
||||
|
||||
%% Compile GIF
|
||||
g.compile(fig.Number);
|
||||
|
||||
disp(fullfile(g.OutputDir, sprintf('%s_fig_%d.gif', g.Name, fig.Number)));
|
||||
@@ -1,55 +1,77 @@
|
||||
% Gitter für lambda0 und S0
|
||||
lambda0_vec = linspace(1260,1360,200);
|
||||
S0_vec = linspace(0.06,0.1,200);
|
||||
[Lambda0, S0] = meshgrid(lambda0_vec, S0_vec);
|
||||
|
||||
% Festen Betriebsparameter
|
||||
lambda = 1293; % nm
|
||||
lambda = 1290; % nm
|
||||
L = 1; % km
|
||||
|
||||
% Dispersion berechnen (lineare Näherung)
|
||||
D = S0 .* ( lambda - Lambda0 ) * L;
|
||||
mu_zwd = 1317; % nm
|
||||
sigma_zwd = 2; % nm
|
||||
mu_s0 = 0.0872; % ps / nm2 km
|
||||
sigma_s0 = 0.0012; % ps / nm2 km
|
||||
rho = -0.5; % Korrelation
|
||||
|
||||
% Gitter für lambda0 und S0
|
||||
lambda0_vec = linspace(mu_zwd-10, mu_zwd+10, 100);
|
||||
S0_vec = linspace(mu_s0-0.01, mu_s0+0.01, 100);
|
||||
|
||||
% Korrigierte meshgrid Reihenfolge
|
||||
[S0, Lambda0] = meshgrid(S0_vec, lambda0_vec);
|
||||
|
||||
% Dispersion berechnen
|
||||
D = (S0./4) .* ( lambda - (Lambda0.^4)./(lambda^3) ) * L;
|
||||
|
||||
%% 2D-Konturplot nur mit Linien und Text
|
||||
%% 2D-Konturplot
|
||||
figure('Color','w');
|
||||
hold on
|
||||
|
||||
% Konturlinien
|
||||
numLevels = 10;
|
||||
levels = linspace(min(D(:)), max(D(:)), numLevels);
|
||||
[C,h] = contour(S0, Lambda0, D, levels, ...
|
||||
'LineWidth',1.5, ...
|
||||
'ShowText','on', ...
|
||||
'LabelFormat','%0.1f');
|
||||
% --- 1. Hintergrund: Dispersions-Konturlinien (Gerundet für TikZ) ---
|
||||
numLevels_D = 10;
|
||||
% Erzeuge glatte, auf 1 Nachkommastelle gerundete Werte
|
||||
|
||||
% cbrewer2-Colormap für die Linien
|
||||
cmap = cbrewer2('div','RdYlGn', numLevels);
|
||||
colormap(cmap);
|
||||
levels_D = round(linspace(min(D(:)), max(D(:)), numLevels_D), 1);
|
||||
levels_D = unique(levels_D); % Falls durch Rundung doppelte Werte entstehen
|
||||
|
||||
% Achsenlinien
|
||||
% yline(1310, '--k','ZDW_{mean}','LabelVerticalAlignment','top','LabelHorizontalAlignment','center');
|
||||
x0 = 0.09;
|
||||
% xline(x0, '--k','S_{0}','LabelHorizontalAlignment','left');
|
||||
% Colormap in der Länge der verbliebenen Level erstellen
|
||||
cmap_bg = cbrewer2('Blues', length(levels_D));
|
||||
|
||||
for i = 1:length(levels_D)
|
||||
% WICHTIG: Das Level als [Wert, Wert] übergeben!
|
||||
[C,h] = contour(S0, Lambda0, D, [levels_D(i), levels_D(i)], ...
|
||||
'Color', cmap_bg(i,:),...
|
||||
'LineWidth', 1.5, ...
|
||||
'ShowText', 'on', ...
|
||||
'LabelFormat', '%0.1f');
|
||||
h.LabelColor = [0,0,0];
|
||||
end
|
||||
|
||||
% --- NEU: Parameter für die Verteilungen ---
|
||||
mu_zwd = 1317; % nm
|
||||
sigma_zwd = 2; % nm
|
||||
mu_s0 = 0.0872; % ps / nm2 km
|
||||
sigma_s0 = 0.0012; % ps / nm2 km
|
||||
rho = -0.5; % Korrelation
|
||||
|
||||
% --- 3. Randverteilungen (1D Gauss) an den Achsen ---
|
||||
s0_vals = linspace(min(S0_vec), max(S0_vec), 500);
|
||||
gauss_s0 = exp(-0.5*((s0_vals - mu_s0)/sigma_s0).^2);
|
||||
scale_s0 = 4; % Skalierung für die Höhe in der Ansicht
|
||||
plot(s0_vals, min(lambda0_vec) + gauss_s0 * scale_s0, 'k', 'LineWidth', 2);
|
||||
|
||||
% Gaussian auf der x-Linie (S0 = 0.09)
|
||||
mu_zwd = 1310; % nm
|
||||
sigma_zwd = 2; % nm
|
||||
zwd_vals = linspace(min(lambda0_vec), max(lambda0_vec), 500);
|
||||
% PDF berechnen
|
||||
gauss_pdf = (1/(sigma_zwd*sqrt(2*pi))) * exp(-0.5*((zwd_vals-mu_zwd)/sigma_zwd).^2);
|
||||
% Normieren und auf eine sichtbare Breite skalieren
|
||||
scale = 0.005; % passt die Maximal-Auslenkung in x-Richtung an
|
||||
x_gauss = x0 + (gauss_pdf/max(gauss_pdf)) * scale;
|
||||
gauss_zwd = exp(-0.5*((zwd_vals - mu_zwd)/sigma_zwd).^2);
|
||||
scale_zwd = 0.005; % Skalierung für die Auslenkung in der Ansicht
|
||||
plot(min(S0_vec) + gauss_zwd * scale_zwd, zwd_vals, 'k', 'LineWidth', 2);
|
||||
|
||||
% Plot
|
||||
% plot(x_gauss, zwd_vals, 'LineWidth',2);
|
||||
% Hilfslinien für die Mittelwerte
|
||||
xline(mu_s0, '--k', 'Alpha', 0.4);
|
||||
yline(mu_zwd, '--k', 'Alpha', 0.4);
|
||||
|
||||
% Achsenbeschriftung & Titel
|
||||
% Achsenbeschriftung & Titel
|
||||
xlabel('S0 [ps / nm2 km]', 'FontSize', 12);
|
||||
% --- Achsenbeschriftung, Titel & Formatierung ---
|
||||
xlabel('S0 ', 'FontSize', 12);
|
||||
ylabel('ZDW [nm]', 'FontSize', 12);
|
||||
title (sprintf('Dispersion: %d km; %d nm', L, lambda), 'FontSize', 14);
|
||||
% title(sprintf('Dispersion: %d km; %d nm', L, lambda), 'FontSize', 14);
|
||||
|
||||
axis([min(S0_vec) max(S0_vec) min(lambda0_vec) max(lambda0_vec)]);
|
||||
grid on
|
||||
hold off
|
||||
hold off
|
||||
|
||||
%% Für den LaTeX Export
|
||||
mat2tikz_improved("C:/Users/Silas/Documents/6971e0b65b380ca6d71c837f/02_IMDD_System/tikz/dispersion/dispersion_contour.tikz")
|
||||
@@ -1,146 +0,0 @@
|
||||
%% ------------------------------------------------------------
|
||||
% Contour plot: λ_null as function of bandwidth (f_target) and reach (L)
|
||||
% ------------------------------------------------------------
|
||||
|
||||
% Parameters
|
||||
lambda0 = 1310e-9; % [m]
|
||||
S0 = 0.08; % [ps/(nm²·km)]
|
||||
c = physconst('lightspeed');
|
||||
|
||||
% Sweep dimensions
|
||||
f_targets = linspace(50e9, 120e9, 100); % [Hz] (x-axis)
|
||||
L_values = linspace(0.5e3, 10e3, 100); % [m] (y-axis)
|
||||
|
||||
lambda_surface = zeros(numel(L_values), numel(f_targets));
|
||||
Dacc_surface = zeros(numel(L_values), numel(f_targets));
|
||||
|
||||
% Outer loop over fiber length (since L must be scalar)
|
||||
for iL = 1:numel(L_values)
|
||||
L = L_values(iL);
|
||||
[lambda_vec, Dacc_vec] = lambda_for_first_null_full(f_targets, L, lambda0, S0);
|
||||
lambda_vec = 2*abs(lambda0 - lambda_vec);
|
||||
|
||||
if 0
|
||||
fprintf('\n- %d km ------------------------------------\n',L);
|
||||
fprintf(' f_null [GHz] lambda [nm] Dacc [ps/nm]\n');
|
||||
fprintf('----------------------------------------------\n');
|
||||
fprintf('%10.1f %8.2f %+8.3f\n',[f_targets(:)/1e9, lambda_vec(:)*1e9, Dacc_vec(:)].');
|
||||
fprintf('----------------------------------------------\n\n');
|
||||
end
|
||||
|
||||
lambda_surface(iL, :) = lambda_vec; % λ for each f_target
|
||||
Dacc_surface(iL, :) = Dacc_vec; % corresponding accumulated dispersion
|
||||
end
|
||||
|
||||
% Convert for plotting
|
||||
lambda_surface_nm = lambda_surface * 1e9; % [nm]
|
||||
L_km = L_values / 1000; % [km]
|
||||
f_GHz = f_targets / 1e9; % [GHz]
|
||||
|
||||
%% Contour plot
|
||||
figure('Color','w');
|
||||
|
||||
% Define wavelength contour levels [nm]
|
||||
lambda_levels = [1260:10:1290, 1290:5:1300, 1300:2.5:1310];
|
||||
lambda_levels = [100:-20:50, 50:-10:30,30:-5:0];
|
||||
|
||||
% Contour plot
|
||||
contour(f_GHz, L_km, lambda_surface_nm, lambda_levels, ...
|
||||
'LineWidth', 1.5, ...
|
||||
'ShowText', 'on', ...
|
||||
'LabelFormat', '%.0f nm');
|
||||
|
||||
% Colormap and colorbar
|
||||
colormap((cbrewer2('RdYlGn',100)));
|
||||
colorbar;
|
||||
clim([0 100]);
|
||||
|
||||
% Axis formatting
|
||||
xlabel('Signal Bandwidth [GHz]');
|
||||
ylabel('Fiber length [km]');
|
||||
legend('$\Delta \lambda$')
|
||||
|
||||
% X-axis ticks at 56 : 16 : 150 GHz
|
||||
xticks(56:8:150);
|
||||
|
||||
grid on; box on;
|
||||
|
||||
|
||||
%% Optional: overlay accumulated-dispersion contours
|
||||
if 0
|
||||
hold on;
|
||||
[CS, h] = contour(f_GHz, L_km, Dacc_surface, 10, 'k--', 'LineWidth', 0.8);
|
||||
clabel(CS, h, 'Color','k', 'FontSize',8);
|
||||
end
|
||||
|
||||
function [lambda_vec, Dacc_vec] = lambda_for_first_null_full(f_target, L, lambda0, S0)
|
||||
% lambda_for_first_null_full (stable, single-branch + validity checks)
|
||||
% --------------------------------------------------------------------
|
||||
% Computes the wavelength(s) at which the first IM/DD fading null
|
||||
% occurs at frequency/ies f_target using the full dispersion model:
|
||||
%
|
||||
% D(lambda) = (S0/4)*(lambda - lambda0^4 / lambda^3)
|
||||
%
|
||||
% Restricted to the NORMAL-dispersion branch (λ < λ0),
|
||||
% and valid only in the O-band (1260–1360 nm).
|
||||
%
|
||||
% Inputs:
|
||||
% f_target - scalar or vector of target null frequencies [Hz]
|
||||
% L - fiber length [m]
|
||||
% lambda0 - zero-dispersion wavelength (ZDW) [m]
|
||||
% S0 - dispersion slope at ZDW [ps/(nm²·km)]
|
||||
%
|
||||
% Outputs:
|
||||
% lambda_vec - wavelength(s) [m] where first null occurs (clamped to O-band)
|
||||
% Dacc_vec - accumulated dispersion(s) [ps/nm] (NaN if out of valid range)
|
||||
% --------------------------------------------------------------------
|
||||
|
||||
c = physconst('lightspeed');
|
||||
S0_si = S0 * 1e3; % ps/(nm²·km) -> s/(m³)
|
||||
|
||||
% Define O-band boundaries (in meters)
|
||||
lambda_min = 1255e-9;
|
||||
lambda_max = 1361e-9;
|
||||
|
||||
% Force column vector
|
||||
f_target = f_target(:);
|
||||
N = numel(f_target);
|
||||
|
||||
lambda_vec = NaN(N,1);
|
||||
Dacc_vec = NaN(N,1);
|
||||
|
||||
for k = 1:N
|
||||
RHS = c * 0.5 / (f_target(k)^2 * L);
|
||||
|
||||
% Normal-dispersion branch (λ < λ0)
|
||||
fun = @(lambda) -(S0_si/4).*(lambda - (lambda0^4)./(lambda.^3)).*lambda.^2 - RHS;
|
||||
|
||||
% Limit the search to [λ_min, λ0)
|
||||
try
|
||||
lambda_sol = fzero(fun, [lambda_min, lambda0 * 0.999]);
|
||||
catch
|
||||
% If the zero is not within bounds, skip this point
|
||||
lambda_sol = NaN;
|
||||
end
|
||||
|
||||
% Validate solution
|
||||
if isnan(lambda_sol) || lambda_sol < lambda_min || lambda_sol > lambda_max
|
||||
lambda_vec(k) = NaN;
|
||||
Dacc_vec(k) = NaN;
|
||||
continue
|
||||
end
|
||||
|
||||
% Compute D(lambda) and accumulated dispersion
|
||||
D_lambda = (S0_si/4) * (lambda_sol - (lambda0^4)/(lambda_sol^3)) / 1e-6; % ps/(nm·km)
|
||||
Dacc_val = D_lambda * (L/1000); % ps/nm
|
||||
|
||||
% Sanity bound on dispersion (avoid unphysical > ±100 ps/nm)
|
||||
if abs(Dacc_val) > 100
|
||||
lambda_vec(k) = NaN;
|
||||
Dacc_vec(k) = NaN;
|
||||
else
|
||||
lambda_vec(k) = lambda_sol;
|
||||
Dacc_vec(k) = Dacc_val;
|
||||
end
|
||||
end
|
||||
end
|
||||
@@ -5,19 +5,34 @@ arguments
|
||||
end
|
||||
cleanfigure;
|
||||
matlab2tikz(char(filename), ...
|
||||
'width','\fwidth', ...
|
||||
'height','\fheight', ...
|
||||
'showInfo',false, ...
|
||||
'extraAxisOptions',{ ...
|
||||
'legend style={font=\footnotesize}', ...
|
||||
'xlabel style={font=\color{white!15!black},font=\small},',...
|
||||
'ylabel style={font=\color{white!15!black},font=\small},',...
|
||||
'width', '\fwidth', ...
|
||||
'height', '\fheight', ...
|
||||
'showInfo', false, ...
|
||||
'extraAxisOptions', { ...
|
||||
'xlabel style={font=\color{white!15!black}\small}', ...
|
||||
'ylabel style={font=\color{white!15!black}\small}', ...
|
||||
'scaled ticks=false', ...
|
||||
'tick label style={/pgf/number format/fixed, /pgf/number format/1000 sep={}}', ...
|
||||
'every axis/.append style={font=\scriptsize}', ...
|
||||
'legend columns=1', ...
|
||||
'every axis/.append style={font=\scriptsize}',...
|
||||
'legend columns=1',...
|
||||
'legend style={at={(0.02,0.98)},font=\footnotesize,draw=black!60,rounded corners=2pt,inner sep=1pt,fill=white,column sep=6pt,anchor= north west}',...
|
||||
'legend style={at={(0.02,0.98)},draw=white!0!white,font=\scriptsize,inner sep=0.1pt,fill=white,column sep=1pt,anchor= north west}',...
|
||||
'every axis/.append style={font=\scriptsize}',...
|
||||
'legend style={at={(0.02,0.98)}, anchor=north west, font=\scriptsize, draw=black!100, rounded corners=2pt, inner sep=1pt, fill=white, column sep=2pt}' ...
|
||||
});
|
||||
% matlab2tikz(char(filename), ...
|
||||
% 'width','\fwidth', ...
|
||||
% 'height','\fheight', ...
|
||||
% 'showInfo',false, ...
|
||||
% 'extraAxisOptions',{ ...
|
||||
% 'legend style={font=\footnotesize}', ...
|
||||
% 'xlabel style={font=\color{white!15!black},font=\small},',...
|
||||
% 'ylabel style={font=\color{white!15!black},font=\small},',...
|
||||
% 'legend columns=1', ...
|
||||
% 'every axis/.append style={font=\scriptsize}',...
|
||||
% 'legend columns=1',...
|
||||
% 'legend style={at={(0.02,0.98)},font=\footnotesize,draw=black!60,rounded corners=2pt,inner sep=1pt,fill=white,column sep=6pt,anchor= north west}',...
|
||||
% 'legend style={at={(0.02,0.98)},draw=white!0!white,font=\scriptsize,inner sep=0.1pt,fill=white,column sep=1pt,anchor= north west}',...
|
||||
% 'every axis/.append style={font=\scriptsize}',...
|
||||
% 'scaled ticks=false', ...
|
||||
% 'tick label style={/pgf/number format/fixed, /pgf/number format/1000 sep={}}' ...
|
||||
% });
|
||||
|
||||
end
|
||||
Reference in New Issue
Block a user