silas/imdd_silas
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge branch 'main' of cau-git.rz.uni-kiel.de:nt/mitarbeiter/silas/imdd_simulation

Browse Source
This commit is contained in:
Silas Oettinghaus
2026-03-25 09:33:44 +01:00
parent e16806b520 5b90bc11a0
commit 69578b3c43
14 changed files with 1033 additions and 24 deletions
Download Patch File Download Diff File Expand all files Collapse all files

6
Classes/00_signals/Signal.m
View File

@@ -925,9 +925,10 @@ classdef Signal
M
options.fignum = 100;
options.displayname = "";
options.mode = 1; %1= histogram method; 2= intuitive "line based" eye
end
mode = 2;
mode = options.mode;
histpoints = 2048; %% verticale resolution
histpoints = floor(histpoints/2)*2+1; %% to have the eye digram centered around one point make the vertical resolution uneven
@@ -966,8 +967,9 @@ classdef Signal
col = cbrewer2('Set1',2);
for n=1:1000
hold on
plot(eye_mat(:,n),'LineStyle',':','LineWidth',0.1,'Color',col(2,:));
plot(eye_mat(:,n),'LineStyle','-','LineWidth',0.1,'Color',col(2,:));
end
xlabel('Samples','Interpreter','latex')
ylabel('Amplitude of Signal','Interpreter','latex');
xlim([0 histpoints_horizontal])

12
Functions/Minimal_examples/matched_filter_minimal.m
View File

@@ -52,11 +52,17 @@ symbols = PAMmapper(M,0).map(bits);
symbols.fs = fsym;
symbols.spectrum("displayname",'Symbols','fignum',1);
%% RRC Shaping
rcalpha = 1;
Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"alpha",rcalpha);
for rcalpha = 0.1:0.2:1
% rcalpha = 0.5;
Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rc","pulselength",16,"alpha",rcalpha);
Digi_sig = Pform.process(symbols);
Digi_sig.spectrum("displayname",'Signal after pluse shaping','fignum',1);
% Digi_sig.spectrum("displayname",'Signal after pluse shaping','fignum',1);
Digi_sig.eye(fsym,M,"fignum",0.1*10,"mode",1);
end
%% RRC Matched Filtering
Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"alpha",rcalpha);

8
Functions/Theory/Dissertation/dispersion_around_zdw.m
View File

@@ -49,6 +49,14 @@ for k = 1:size(scenarios, 1)
% 'HandleVisibility', 'off'); % Hide from legend
% Capture Wide Spec (k=1) bounds for the TikZ measurement lines
hp.FaceAlpha = 0.5;
% 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

116
Functions/Theory/Dissertation/pmd_vs_length.m Normal file
View File

@@ -0,0 +1,116 @@
% pmd_vs_length.m
% ------------------------------------------------------------
% Plots the Foschini-Poole (1991) analytical variance formula for PMD:
%
% sigma_T^2(z) = 2*(Delta_beta1)^2 * lc^2
% * [ exp(-z/lc) + z/lc - 1 ]
%
% and overlays the two asymptotic regimes:
% - Short-reach (z << lc) : sigma_T(z) ~ (Delta_beta1) * z
% - Long-haul (z >> lc) : sigma_T(z) ~ Dp * sqrt(z)
%
% Parameters follow typical SMF values from the literature.
% ------------------------------------------------------------
clear; clc;
%% ── Parameters ──────────────────────────────────────────────────────────────
% Intrinsic local birefringence [ps/km]
Delta_beta1 = 1e-1; % typical value, adjust as needed
% Correlation length [km]
lc = 0.01; % ~50 m, typical for G.652 SMF
% PMD parameter [ps / sqrt(km)] — derived from the two above
Dp = Delta_beta1 * sqrt(2 * lc);
% Distance axis [km]
z_max = 10; % maximum distance
z = linspace(0.001, z_max, 10000); % avoid z = 0 in log plot
%% ── Exact Foschini-Poole formula (sigma_T in ps) ────────────────────────────
sigma_T_sq = 2 .* Delta_beta1.^2 .* lc.^2 ...
.* (exp(-z ./ lc) + z ./ lc - 1);
sigma_T = sqrt(sigma_T_sq); % RMS DGD [ps]
%% ── Asymptotic regimes ───────────────────────────────────────────────────────
% Short-reach: linear growth (z << lc)
sigma_T_short = Delta_beta1 .* z; % [ps]
% Long-haul: square-root growth (z >> lc)
sigma_T_long = Dp .* sqrt(z); % [ps]
%% ── Plot ─────────────────────────────────────────────────────────────────────
figure('Color','w','Position',[100 100 760 480]);
hold on;
% Color palette (matching dissertation style)
c_exact = [0.1216, 0.4706, 0.7059]; % blue – exact
c_short = [0.8392, 0.1529, 0.1569]; % red – short-reach asymptote
c_long = [0.1961, 0.6314, 0.1725]; % green – long-haul asymptote
% Exact solution
h_exact = plot(z, sigma_T, ...
'Color', c_exact, 'LineWidth', 2.0, ...
'DisplayName', 'Exact (Foschini \& Poole)');
% Short-reach asymptote σ_T ≈ Δβ₁ · z
h_short = plot(z, sigma_T_short, ...
'Color', c_short, 'LineWidth', 1.4, 'LineStyle', '--', ...
'DisplayName', '$\sigma_T \approx \Delta\beta_1 \cdot z$ \quad ($z \ll l_c$)');
% Long-haul asymptote σ_T ≈ D_p √z
h_long = plot(z, sigma_T_long, ...
'Color', c_long, 'LineWidth', 1.4, 'LineStyle', ':', ...
'DisplayName', '$\sigma_T \approx D_p \sqrt{z}$ \quad ($z \gg l_c$)');
%% ── Axes & decoration ────────────────────────────────────────────────────────
ax = gca;
set(ax, 'XScale', 'log', 'YScale', 'log');
% ── X-axis: linear-style tick labels on log scale ────────────────────────
x_ticks = [1e-3, 1e-2, 1e-1, 1, 10];
ax.XTick = x_ticks;
ax.XTickLabel = arrayfun(@(v) sprintf('%g km', v), x_ticks, 'UniformOutput', false);
% ── Y-axis: linear-style tick labels on log scale ────────────────────────
y_ticks = [1e-3, 1e-2, 1e-1, 1, 10];
ax.YTick = y_ticks;
ax.YTickLabel = arrayfun(@(v) sprintf('%g ps', v), y_ticks, 'UniformOutput', false);
xlabel('Fiber length $z$ [km]', 'Interpreter', 'latex');
ylabel('RMS DGD $\sigma_T$ [ps]', 'Interpreter', 'latex');
grid on; box on;
xlim([min(z) z_max]);
legend([h_exact, h_short, h_long], ...
'Location', 'northwest', 'Interpreter', 'latex', 'FontSize', 9);
% Parameter annotation
anno_str = sprintf( ...
['$\\Delta\\beta_1 = %.3g$ ps/km\n' ...
'$l_c = %.0f$ m\n' ...
'$D_p = \\Delta\\beta_1\\sqrt{2l_c} = %.4g$ ps/$\\sqrt{\\mathrm{km}}$'], ...
Delta_beta1, lc*1e3, Dp);
annotation('textbox', [0.57 0.14 0.38 0.22], ...
'String', anno_str, ...
'Interpreter', 'latex', ...
'FontSize', 8.5, ...
'BackgroundColor','w', ...
'EdgeColor', [0.5 0.5 0.5], ...
'LineWidth', 0.8, ...
'FitBoxToText', 'on');
%% ── Regime transition marker ─────────────────────────────────────────────────
% Mark the crossover region around z = lc
xline(lc, '--', ...
'Color', [0.5 0.5 0.5], 'LineWidth', 0.8, ...
'HandleVisibility', 'off');
text(lc * 1.15, min(sigma_T)*3, '$l_c$', ...
'Interpreter', 'latex', 'Color', [0.4 0.4 0.4], 'FontSize', 9);
%% ── Export (uncomment to use) ────────────────────────────────────────────────
% mat2tikz_improved('C:\...\tikz\pmd\pmd_vs_length.tikz')

25
Libs/mutual information rate TUM/LICENSE.txt Normal file
View File

@@ -0,0 +1,25 @@
Copyright (c) 2019 Francisco Javier Garcia-Gomez <javier.garcia@tum.de>
Institute for Communications Engineering (LNT)
Technical University of Munich, Germany
www.lnt.ei.tum.de
All rights reserved.
Permission is hereby granted, free of charge, to any person obtaining a
copy of this software and associated documentation files (the
"Software"), to deal in the Software without restriction, including
without limitation the rights to use, copy, modify, merge, publish,
distribute, sublicense, and/or sell copies of the Software, and to permit
persons to whom the Software is furnished to do so, subject to the
following conditions:
The above copyright notice and this permission notice shall be included
in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN
NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE
USE OR OTHER DEALINGS IN THE SOFTWARE.

24
Libs/mutual information rate TUM/README.md Normal file
View File

@@ -0,0 +1,24 @@
## MI-CG: Numerically Computing Achievable Rates of Memoryless Channels
This repository provides a MATLAB function mi_cg.m to numerically compute achievable rates for memoryless channels. The function uses a conditionally-Gaussian (CG) channel model to obtain a lower bound on the achievable rate of the true channel. The method is well-known, and it is explained in [this short document](https://mediatum.ub.tum.de/node?id=1533663). Two example scripts that compute several achievable rate curves are also provided.
### Citation
This software and the accompanying document are meant as a tutorial to get started with mutual information as a numerical figure of merit for a communications channel. The software is provided under the open-source [MIT license](https://opensource.org/licenses/MIT). If you use the software in your academic work, please cite the accompanying [document](https://mediatum.ub.tum.de/node?id=1533663) as follows:
> F. J. Garcia-Gomez, “Numerically computing achievable rates of memoryless channels,” TUM University Library, 2019. [Online]. Available: https://mediatum.ub.tum.de/node?id=1533663
The corresponding BibTeX entry is
```
@article{garcia2019numerically,
author = "Francisco Javier Garcia-Gomez",
title = "Numerically Computing Achievable Rates of Memoryless Channels",
year = "2019",
journal="TUM University Library",
url={https://mediatum.ub.tum.de/node?id=1533663}
}
```
### Acknowledgment
This work was supported by the German Research Foundation (DFG) under Grant KR 3517/8-2.

156
Libs/mutual information rate TUM/air.m Normal file
View File

@@ -0,0 +1,156 @@
function air = air(x,r,idx_tx,Px,M_training)
MAX_MEMORY = 200e6; % maximum allowed size for a matrix
if nargin == 3
Px = [];
M_training = [];
end
if nargin == 4
M_training = [];
end
% if input is complex, separate into real and imaginary parts
if any(imag(x(:))~=0) || any(imag(r(:))~=0)
x = [real(x); imag(x)];
r = [real(r); imag(r)];
end
D = size(x, 1); % D = 2 if complex x
N = size(x, 2); % number of constellation points
M = size(r, 2); % number of samples
% set default training set size
if isempty(M_training)
M_training = ceil(0.3*M);
end
M_testing = M - M_training;
% Training: estimate parameters of the conditionally Gaussian model
% sort according to transmit index
[idx_tx_training, idx_sort] = sort(idx_tx(1:M_training));
r_training = r(:, idx_sort);
i_bounds = zeros(1, N+1);
% compute conditional means and covariance matrices
C_n = zeros(D, D, N);
det_n = zeros(1, N);
for n=1:N
% find how many times x(:, n) was transmitted and update i_bounds
N_current_x = find(idx_tx_training((i_bounds(n)+1):end)==n, 1, 'last');
if isempty(N_current_x), N_current_x=0; end
i_bounds(n+1) = i_bounds(n) + N_current_x;
if N_current_x > 0
% Compute mu_n=E[Y|X=x_n] according to Eq. (14) and store it in
% x(:, n) to save space
x(:, n) = sum(r_training(:, (i_bounds(n)+1):i_bounds(n+1)), 2)/(i_bounds(n+1)-i_bounds(n));
% compute C_n=cov[Y|X=x_n] according to Eq. (15)
r_meanfree = r_training(:, (i_bounds(n)+1):i_bounds(n+1)) - x(:, n);
C_n(:, :, n) = (r_meanfree*r_meanfree')/(i_bounds(n+1)-i_bounds(n));
% store also the determinant of C(:, :, n)
det_n(n) = det(C_n(:, :, n));
% if the determinant is 0, or if the matrix is badly conditioned,
% regularize by adding a small identity matrix. Note that we do
% need the check for 0 determinant, in case a cloud has exactly 0
% variance according to the training set
if det_n(n)==0 || cond(C_n(:, :, n))>1e16
C_n(:, :, n) = C_n(:, :, n) + 5 * eps * eye(D);
det_n(n) = (5*eps)^D;
end
end
end
% uniform input pmf Px if not provided
if isempty(Px)
Px = repmat(1/N, [1, N]);
end
% extract testing set and sort it according to transmit index
[idx_tx_testing, idx_sort] = sort(idx_tx((M_training+1):M));
r_testing = r(:, M_training+idx_sort);
% computation of h(Y|X)
h_Y_X = 0;
i_bounds_testing = zeros(1, N+1);
% loop over constellation points to compute h(Y|X)
for n = 1:N
% find how many times x(:, n) was transmitted and update
% i_bounds_testing
N_current_x = find(idx_tx_testing((i_bounds_testing(n)+1):end)==n, 1, 'last');
if isempty(N_current_x), N_current_x=0; end
i_bounds_testing(n+1) = i_bounds_testing(n) + N_current_x;
% add the corresponding contribution to the mutual information (two
% first lines of Eq. (17)). This, together with
% D/2*log2(2*pi) after the end of the loop, gives h(Y|X)
h_Y_X = h_Y_X + N_current_x * log2(det_n(n))/2+...
sum(sum(conj(r_testing(:, (i_bounds_testing(n)+1):i_bounds_testing(n+1))-x(:, n)).*(C_n(:, :, n)\(r_testing(:, (i_bounds_testing(n)+1):i_bounds_testing(n+1))-x(:, n)))))/2/log(2);
end
h_Y_X = D/2*log2(2*pi) + h_Y_X/M_testing;
% When computing log(py), we might run out of memory. If necessary, we
% doe the computation in blocks
logpy = zeros(1, M_testing);
BLOCK_SIZE = floor(MAX_MEMORY/N);
N_blocks = ceil(M_testing/BLOCK_SIZE);
% loop over blocks of symbols. This loop can be replaced by parfor to allow
% parallel computation
for i_block = 1:N_blocks
logpy_cur = zeros(1, M_testing);
% beginning of block
i_start = (i_block-1) * BLOCK_SIZE + 1;
% end of block
i_end = min(M_testing, i_block*BLOCK_SIZE);
% block size
current_block_size = i_end-i_start+1;
% compute exponents of third line of (17)
exponents = zeros(N, current_block_size);
for n = 1:N
exponents(n, :) = -log(det_n(n))/2-real(sum(conj(r_testing(:, i_start:i_end)-x(:, n)).*(C_n(:, :, n)\(r_testing(:, i_start:i_end)-x(:, n))), 1))/2;
%sum über 2 einträge von r
end
% compute third line of Eq. (17). Use a custom function
% that computes log(sum(exp(x))) avoiding overflow errors
logpy_cur(i_start:i_end) = math_logsumexp(log(Px(:))+exponents, 1);
logpy = logpy + logpy_cur;
end
% output entropy h(Y)
h_Y = D/2*log2(2*pi) - mean(logpy)/log(2);%log basis change
% compute mutual information
air = h_Y - h_Y_X;
end
function [y] = math_logsumexp(x, dim)
%[y] = math_logsumexp(x, dim)
% Computes log(sum(exp(x), dim)), avoiding overflow errors when one of the
% x is large.
if nargin<2 || isempty(dim)
m = max(x);
y = m + log(sum(exp(x-m)));
else
m = max(x, [], dim);
y = m + log(sum(exp(x-m), dim));
end
end

124
Libs/mutual information rate TUM/example_mi_cg_complex_16qam.m Normal file
View File

@@ -0,0 +1,124 @@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright (c) 2019 Francisco Javier Garcia-Gomez <javier.garcia@tum.de>
% Institute for Communications Engineering (LNT)
% Technical University of Munich, Germany
% www.lnt.ei.tum.de
%
% All rights reserved.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Permission is hereby granted, free of charge, to any person obtaining a
% copy of this software and associated documentation files (the
% "Software"), to deal in the Software without restriction, including
% without limitation the rights to use, copy, modify, merge, publish,
% distribute, sublicense, and/or sell copies of the Software, and to permit
% persons to whom the Software is furnished to do so, subject to the
% following conditions:
%
% The above copyright notice and this permission notice shall be included
% in all copies or substantial portions of the Software.
%
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
% OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
% MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN
% NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
% DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
% OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE
% USE OR OTHER DEALINGS IN THE SOFTWARE.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%% example_mi_cg_complex_16qam.m %%%%%%%%%%%%%%%%%%%%%%
%
% Example usage of the function mi_cg to numerically compute mutual
% information between two complex sequences. This file simulates a
% one-dimensional complex AWGN channel with 16-QAM constellation for
% different SNRs, and then plots the achievable rate, which is equal to two
% times the 4-PAM curve of Fig. 1 of [1].
%
% Note that using mi_cg with complex sequences assumes that the channel
% model q(Y|X) is circularly symmetric for a given X=x (i.e., that the
% received clouds are circular). This is not the case in a channel with
% phase noise: see example_it_mi_cg.m for an example of how to deal with
% non-circularly-symmetric channels.
%
% [1] G. David Forney and Gottfried Ungerboeck, "Modulation and Coding for
% Linear Gaussian Channels", IEEE Trans. Inf. Theory vol. 44, no. 6, pp.
% 2384-2415, October 1998
%
% Technische Universitaet Muenchen - Lehrstuhl fuer Nachrichtentechnik
% Date: 12.12.2019
% Author: Francisco Javier Garcia-Gomez <javier.garcia@tum.de>
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
% close all;
%%% Simulation parameters
M_QAM=16; % QAM size
var_w=1; % noise variance
SNR_dB_values=(-5):2:25; % SNR values in dB
N=8000; % number of Monte-Carlo points
%%% Derived parameters
n_SNR=length(SNR_dB_values);
SNR_values=10.^(SNR_dB_values/10);
%%% Generation of QAM constellation with power 1
M_PAM=sqrt(M_QAM); % number of points per real dimension
X_1d=-(M_PAM-1)+2*(0:(M_PAM-1)); % generate PAM
X=X_1d+1i*X_1d.'; % transform to QAM
X=X(:).'*sqrt(3/2/(M_QAM-1)); % set to unit power
X=[real(X); imag(X)]; % separate real and imaginary parts
%%% Uniformly choose transmit indices
idx_tx=randi(M_QAM, [1, N]);
%%% AWGN noise
w=sqrt(var_w)*(randn([2, N]));%+1i*randn([1, N]));
%%% Loop over the SNR values
MI_awgn=zeros(1, n_SNR);
fig = figure(1);
for i_SNR=1:n_SNR
% transmitted points
s=sqrt(SNR_values(i_SNR)*var_w)*X(:, idx_tx);
% AWGN channel
r_awgn = s+w;
%%%%%%%%%%%%
clf
hold on
xlim([-20, 20]);
ylim([-20, 20]);
%received signal
scatter(r_awgn(1,:),r_awgn(2,:),1,"red",'.');
%tx signal constellation
scatter(s(1,:),s(2,:),4,"black",'o','filled');
%uni power transmit constallation
scatter(X(1,:),X(2,:),5,'blue','o','filled');
drawnow
pause(0.1)
%%%%%%%%%%%%
% Compute MI
MI_awgn(i_SNR)=air(X, r_awgn, idx_tx);
% The following also works but is slower
% MI_awgn(i_SNR)=mi_cg(s, r_awgn);
end
I_shannon=log2(1+SNR_values);
figure(2);
hold on
plot(SNR_dB_values, I_shannon, '-', 'DisplayName', 'log_2 (1+SNR)');
hold on;
plot(SNR_dB_values, MI_awgn, '--', 'DisplayName', [num2str(M_QAM) '-QAM, AWGN']);
hold off;
xlabel('SNR (dB)'); ylabel('Achievable rate (bits/complex dimension)');
title(['Achievable rate of ' num2str(M_QAM) '-QAM in AWGN']);
legend('Location', 'NorthWest');

140
Libs/mutual information rate TUM/example_mi_cg_real_2D_awgn_vs_phasenoise.m Normal file
View File

@@ -0,0 +1,140 @@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright (c) 2019 Francisco Javier Garcia-Gomez <javier.garcia@tum.de>
% Institute for Communications Engineering (LNT)
% Technical University of Munich, Germany
% www.lnt.ei.tum.de
%
% All rights reserved.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Permission is hereby granted, free of charge, to any person obtaining a
% copy of this software and associated documentation files (the
% "Software"), to deal in the Software without restriction, including
% without limitation the rights to use, copy, modify, merge, publish,
% distribute, sublicense, and/or sell copies of the Software, and to permit
% persons to whom the Software is furnished to do so, subject to the
% following conditions:
%
% The above copyright notice and this permission notice shall be included
% in all copies or substantial portions of the Software.
%
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
% OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
% MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN
% NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
% DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
% OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE
% USE OR OTHER DEALINGS IN THE SOFTWARE.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%% example_mi_cg_real_2D_awgn_vs_phasenoise %%%%%%%%%%%%%%%%%
%
% Example usage of the function mi_cg to numerically compute mutual
% information between two sequences. This file simulates a two-dimensional
% AWGN channel with 4-PAM constellation (or, equivalently, a complex AWGN
% channel with 16-QAM constellation) for different SNRs, and then plots
% the achievable rate, reproducing the 4-PAM curve of Fig. 1 of [1]. Note
% that this curve can also be reproduced by simulating a 1-dimensional
% channel with 4-PAM: this file is an example of how to deal with multiple
% dimensions.
%
% This file also simulates a complex phase-noise channel:
% r=s.*exp(1i*theta)+w
% where w is complex AWGN and theta is i.i.d. real Gaussian. As q(Y|X) for
% this channel is not circularly-symmetric, this complex channel needs to
% be separated into two real dimensions to compute the mutual information.
% The resulting achievable rate is, as expected, below the rate for the
% AWGN channel.
%
% [1] G. David Forney and Gottfried Ungerboeck, "Modulation and Coding for
% Linear Gaussian Channels", IEEE Trans. Inf. Theory vol. 44, no. 6, pp.
% 2384-2415, October 1998
%
% [2]
%
% Technische Universitaet Muenchen - Lehrstuhl fuer Nachrichtentechnik
% Date: 12.12.2019
% Author: Francisco Javier Garcia-Gomez <javier.garcia@tum.de>
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
close all;
%%% Simulation parameters
D=1; % number of complex dimensions
M_QAM=16; % QAM size
var_w=1; % noise variance
E_theta=0.2; % phase offset for the phase noise channel
var_theta=0.05; % phase noise variance
SNR_dB_values=(-5):2:25; % SNR values in dB
N=8000; % number of Monte-Carlo points
%%% Derived parameters
n_SNR=length(SNR_dB_values);
SNR_values=10.^(SNR_dB_values/10);
%%% Generation QAM constellation with power 1
M_PAM=sqrt(M_QAM); % number of points per real dimension
X_1d=-(M_PAM-1)+2*(0:(M_PAM-1)); % generate PAM
X=X_1d+1i*X_1d.'; % transform to QAM
X=X(:).'*sqrt(3/2/(M_QAM-1)); % set to unit power
X=[real(X); imag(X)]; % separate real and imaginary parts
%%% Uniformly choose transmit indices
idx_tx=randi(M_QAM^D, [1, N]);
%%% AWGN noise (real and imaginary parts)
w=sqrt(var_w/2)*randn([2*D, N]);
%%% phase noise
theta=E_theta+sqrt(var_theta)*randn([D, N]);
%%% Loop over the SNR values
MI_awgn=zeros(1, n_SNR);
MI_phasenoise=zeros(1, n_SNR);
fig = figure(1);
for i_SNR=1:n_SNR
% transmitted points
s=sqrt(SNR_values(i_SNR)*var_w)*X(:, idx_tx);
% AWGN channel
r_awgn=s+w;
% Compute MI
MI_awgn(i_SNR)=air(X, r_awgn, idx_tx)/D;
% The following also works but is slower
% MI_awgn(i_SNR)=air(s, r_awgn)/D;
% Phase-noise channel
s_complex=s(1:2:end, :)+1i*s(2:2:end, :); % transform to complex
r_phasenoise_complex=s_complex.*exp(1i*theta); % phase noise and AWGN
r_phasenoise=zeros(2*D, N); % transform to real
r_phasenoise(1:2:end, :)=real(r_phasenoise_complex);
r_phasenoise(2:2:end, :)=imag(r_phasenoise_complex);
r_phasenoise=r_phasenoise+w; % add AWGN
clf
fig = scatter(r_phasenoise(1,:),r_phasenoise(2,:),1,'.');
xlim([-20, 20]);
ylim([-20, 20]);
drawnow
pause(0.1)
% Compute MI
MI_phasenoise(i_SNR)=air(X, r_phasenoise, idx_tx)/D;
% The following also works but is slower
% MI_phasenoise(i_SNR)=air(s, r_phasenoise)/D;
end
% Shannon capacity of the AWGN channel
I_shannon=log2(1+SNR_values);
figure;
plot(SNR_dB_values, I_shannon, '-', 'DisplayName', 'log_2 (1+SNR)');
hold on;
plot(SNR_dB_values, MI_awgn, '--', 'DisplayName', '16-QAM, AWGN');
plot(SNR_dB_values, MI_phasenoise, '-.', 'DisplayName', ['16-QAM, Phase noise, \sigma_\Theta^2=' num2str(var_theta)]);
hold off;
xlabel('SNR (dB)'); ylabel('Achievable rate (bits/complex dimension)');
title(['Achievable rate of 16-QAM in AWGN and in a phase noise channel with \theta~N(' num2str(E_theta) ', ' num2str(var_theta) ')']);
legend('Location', 'NorthWest');

164
Libs/mutual information rate TUM/example_silas.m Normal file
View File

@@ -0,0 +1,164 @@
colored_noise = true;
I_shannon=[];
for cmplx = [0,1]
cnt = 1;
for M_PAM = [4,8,16]
complex_constellation = cmplx;
%%% Simulation parameters
% M_PAM=2; % QAM size
var_w=1; % noise variance
SNR_dB_values=(-5):1:35; % SNR values in dB
N=8000; % number of Monte-Carlo points
%%% Derived parameters
n_SNR=length(SNR_dB_values);
SNR_values=10.^(SNR_dB_values/10);
%%% Generation of QAM constellation with power 1
M=sqrt(M_PAM); % number of points per real dimension
if complex_constellation
X_ = qammod(0:M_PAM-1,M_PAM,"gray");
else
X_ = pammod(0:M_PAM-1,M_PAM,0,'gray');
end
X_ = X_ ./ rms(unique(X_));
X_=[real(X_); imag(X_)];
%%%%%%%%%%%%
figure(10);
clf
hold on
scatter(X_(1,:),X_(2,:),15,'red','x','DisplayName','Matlab');
%%%%%%%%%%%%
%%% Uniformly choose transmit indices
idx_tx=randi(M_PAM, [1, N]);
%%% AWGN noise
w=sqrt(var_w)*(randn([2, N]));
%%% Loop over the SNR values
MI_awgn=zeros(1, n_SNR);
fig = figure(2);
if colored_noise
for i_SNR = 1:n_SNR
% Transmitted signal for the given indices
signal = X_(:, idx_tx);
% Generate white Gaussian noise of the same size as the signal
noise_white = randn(size(signal));
% Define the filter coefficient for colored noise (adjust as needed)
a = 0; % A higher value gives more correlation
% Filter the white noise to create colored noise.
% Here we filter each row (dimension) independently.
noise_colored = zeros(size(signal));
noise_colored(1,:) = filter(1, [1, -a], noise_white(1,:));
noise_colored(2,:) = filter(1, [1, -a], noise_white(2,:));
% Compute signal and unscaled noise power for proper scaling.
signal_power = mean(abs(signal(:)).^2);
noise_power = mean(abs(noise_colored(:)).^2);
% Convert desired SNR from dB to linear scale.
SNR_linear = 10^(SNR_dB_values(i_SNR)/10);
% Scale the colored noise so that signal_power / noise_power equals SNR_linear.
scaling_factor = sqrt(signal_power / (SNR_linear * noise_power));
noise_colored_scaled = scaling_factor * noise_colored;
% Received signal is the sum of the signal and the scaled colored noise.
r_colored = signal + noise_colored_scaled;
% For SNR measurement, compute the noise actually added.
measured_noise = r_colored - signal;
snr_meas_(i_SNR) = snr(signal(1,:) + 1i*signal(2,:), ...
measured_noise(1,:) + 1i*measured_noise(2,:));
% Plotting the transmitted and received signal
clf;
hold on;
xlim([-4, 4]);
ylim([-4, 4]);
scatter(signal(1,:), signal(2,:), 3, "black", 'o', 'DisplayName','Tx Constellation');
scatter(r_colored(1,:), r_colored(2,:), 1, "red", 'x', 'DisplayName', 'Colored Noise Mapping');
drawnow;
% Compute Mutual Information (or any other metric) with the new channel
MI_awgn_(i_SNR) = air(X_, r_colored, idx_tx);
end
else
for i_SNR=1:n_SNR
s_ = sqrt(SNR_values(i_SNR)*var_w) * X_(:, idx_tx);
% r_awgn_ = s_ + w;
r_awgn_ = awgn(X_(:, idx_tx),SNR_dB_values(i_SNR),"measured",10);
w_awgn_matlab = r_awgn_ - X_(:,idx_tx);
snr_meas_(i_SNR) = snr(X_(1, idx_tx)+1i*X_(2, idx_tx) , w_awgn_matlab(1,:)+1i*w_awgn_matlab(2,:));
%%%%%%%%%%%%
clf
hold on
xlim([-4, 4]);
ylim([-4, 4]);
%received signal
scatter(X_(1, idx_tx),X_(2, idx_tx),3,"black",'o','DisplayName','Tx Constellation');
scatter(r_awgn_(1,:),r_awgn_(2,:),1,"red",'x','DisplayName','Matlab Mapping');
drawnow
%%%%%%%%%%%%
MI_awgn_(i_SNR)=air(X_, r_awgn_, idx_tx);
% The following also works but is slower
% MI_awgn(i_SNR)=mi_cg(s, r_awgn);
end
end
figure(3);
hold on
if isempty(I_shannon)
I_shannon=log2(1+db2pow(snr_meas_));
plot(snr_meas_, I_shannon, '-', 'DisplayName', 'log_2 (1+SNR)','Color','black');
end
cols = [ 0.9047 0.1918 0.1988
0.2941 0.5447 0.7494
0.3718 0.7176 0.3612
1.0000 0.5482 0.1000
0.8650 0.8110 0.4330
0.6859 0.4035 0.2412];
if complex_constellation
plot(snr_meas_, MI_awgn_, ':', 'DisplayName', ['',num2str(M_PAM) '-QAM, AWGN'],'LineWidth',1,'Color',cols(cnt,:));
else
plot(snr_meas_, MI_awgn_, '-', 'DisplayName', ['',num2str(M_PAM) '-PAM, AWGN'],'LineWidth',1,'Color',cols(cnt,:));
end
hold off;
ylim([0,8]);
xlim([min(SNR_dB_values) max(SNR_dB_values)])
xlabel('SNR (dB)'); ylabel('Achievable rate (bits/complex dimension)');
title(['Achievable rate of ' num2str(M_PAM) '-QAM in AWGN']);
legend('Location', 'NorthWest');
yline(log2(M_PAM),'LineStyle',':','HandleVisibility','off','Color','black');
cnt = cnt+1;
end
end

104
Libs/mutual information rate TUM/mi_cg.m Normal file
View File

@@ -0,0 +1,104 @@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright (c) 2019 Francisco Javier Garcia-Gomez <javier.garcia@tum.de>
% Institute for Communications Engineering (LNT)
% Technical University of Munich, Germany
% www.lnt.ei.tum.de
%
% All rights reserved.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Permission is hereby granted, free of charge, to any person obtaining a
% copy of this software and associated documentation files (the
% "Software"), to deal in the Software without restriction, including
% without limitation the rights to use, copy, modify, merge, publish,
% distribute, sublicense, and/or sell copies of the Software, and to permit
% persons to whom the Software is furnished to do so, subject to the
% following conditions:
%
% The above copyright notice and this permission notice shall be included
% in all copies or substantial portions of the Software.
%
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
% OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
% MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN
% NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
% DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
% OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE
% USE OR OTHER DEALINGS IN THE SOFTWARE.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%% example_mi_cg_complex_16qam.m %%%%%%%%%%%%%%%%%%%%%%
%
% Example usage of the function mi_cg to numerically compute mutual
% information between two complex sequences. This file simulates a
% one-dimensional complex AWGN channel with 16-QAM constellation for
% different SNRs, and then plots the achievable rate, which is equal to two
% times the 4-PAM curve of Fig. 1 of [1].
%
% Note that using mi_cg with complex sequences assumes that the channel
% model q(Y|X) is circularly symmetric for a given X=x (i.e., that the
% received clouds are circular). This is not the case in a channel with
% phase noise: see example_it_mi_cg.m for an example of how to deal with
% non-circularly-symmetric channels.
%
% [1] G. David Forney and Gottfried Ungerboeck, "Modulation and Coding for
% Linear Gaussian Channels", IEEE Trans. Inf. Theory vol. 44, no. 6, pp.
% 2384-2415, October 1998
%
% Technische Universitaet Muenchen - Lehrstuhl fuer Nachrichtentechnik
% Date: 12.12.2019
% Author: Francisco Javier Garcia-Gomez <javier.garcia@tum.de>
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
close all;
%%% Simulation parameters
M_QAM=16; % QAM size
var_w=1; % noise variance
SNR_dB_values=(-5):2:25; % SNR values in dB
N=8000; % number of Monte-Carlo points
%%% Derived parameters
n_SNR=length(SNR_dB_values);
SNR_values=10.^(SNR_dB_values/10);
%%% Generation of QAM constellation with power 1
M_PAM=sqrt(M_QAM); % number of points per real dimension
X_1d=-(M_PAM-1)+2*(0:(M_PAM-1)); % generate PAM
X=X_1d+1i*X_1d.'; % transform to QAM
X=X(:).'*sqrt(3/2/(M_QAM-1)); % set to unit power
%%% Uniformly choose transmit indices
idx_tx=randi(M_QAM, [1, N]);
%%% AWGN noise
w=sqrt(var_w/2)*(randn([1, N])+1i*randn([1, N]));
%%% Loop over the SNR values
MI_awgn=zeros(1, n_SNR);
for i_SNR=1:n_SNR
% transmitted points
s=sqrt(SNR_values(i_SNR)*var_w)*X(:, idx_tx);
% AWGN channel
r_awgn = s+w;
% Compute MI
MI_awgn(i_SNR)=air(X, r_awgn, idx_tx);
% The following also works but is slower: MI_awgn(i_SNR)=air(s, r_awgn);
end
I_shannon=log2(1+SNR_values);
figure;
plot(SNR_dB_values, I_shannon, '-', 'DisplayName', 'log_2 (1+SNR)');
hold on;
plot(SNR_dB_values, MI_awgn, '--', 'DisplayName', [num2str(M_QAM) '-QAM, AWGN']);
hold off;
xlabel('SNR (dB)'); ylabel('Achievable rate (bits/complex dimension)');
title(['Achievable rate of ' num2str(M_QAM) '-QAM in AWGN']);
legend('Location', 'NorthWest');

135
Libs/mutual information rate TUM/silas_example_mi_cg_pam.m Normal file
View File

@@ -0,0 +1,135 @@
clear;
I_shannon=[];
for cmplx = [1,0]
cnt = 1;
for M_PAM = [2,4,8,16]
% close all;
usegarcia = 0;
complex_constellation = cmplx;
%%% Simulation parameters
% M_PAM=2; % QAM size
var_w=1; % noise variance
SNR_dB_values=(-5):1:35; % SNR values in dB
N=8000; % number of Monte-Carlo points
%%% Derived parameters
n_SNR=length(SNR_dB_values);
SNR_values=10.^(SNR_dB_values/10);
%%% Generation of QAM constellation with power 1
M=sqrt(M_PAM); % number of points per real dimension
if complex_constellation
X_ = qammod(0:M_PAM-1,M_PAM,"gray");
else
X_ = pammod(0:M_PAM-1,M_PAM,0,'gray');
end
X_ = X_ ./ rms(unique(X_));
X_=[real(X_); imag(X_)];
%%%%%%%%%%%%
figure(10);
clf
hold on
scatter(X_(1,:),X_(2,:),15,'red','x','DisplayName','Matlab');
%%%%%%%%%%%%
%%% Uniformly choose transmit indices
idx_tx=randi(M_PAM, [1, N]);
%%% AWGN noise
w=sqrt(var_w)*(randn([2, N]));
%%% Loop over the SNR values
MI_awgn=zeros(1, n_SNR);
fig = figure(2);
for i_SNR=1:n_SNR
if usegarcia
%scale transmitted points acc. to SNR condition (this is not correct I think, the papaer also show diff results)
% s=sqrt(SNR_values(i_SNR)*var_w)*X(:, idx_tx);
%scale nosie acc. to snr condition
noise_power = mean(abs(X_(1, idx_tx)+1i*X_(2, idx_tx)).^2) / SNR_values(i_SNR); % Noise power
w_ = sqrt(noise_power) .* w / sqrt(2);
% AWGN channel
r_awgn = X(:, idx_tx)+w_;
% snr_meas(i_SNR) = snr(s(1,:)+1i*s(2,:),w(1,:)+1i*w(2,:));
snr_meas(i_SNR) = snr(X_(1, idx_tx)+1i*X_(2, idx_tx),w_(1,:)+1i*w_(2,:));
end
s_ = sqrt(SNR_values(i_SNR)*var_w) * X_(:, idx_tx);
% r_awgn_ = s_ + w;
if cmplx
r_awgn_ = awgn(X_(:, idx_tx),SNR_dB_values(i_SNR),"measured",10);
w_awgn_matlab = r_awgn_ - X_(:,idx_tx);
snr_meas_(i_SNR) = snr(X_(1, idx_tx)+1i*X_(2, idx_tx) , w_awgn_matlab(1,:)+1i*w_awgn_matlab(2,:));
else
r_awgn_ = awgn(X_(1, idx_tx),SNR_dB_values(i_SNR),"measured",10);
w_awgn_matlab = r_awgn_ - X_(1,idx_tx);
snr_meas_(i_SNR) = snr(X_(1, idx_tx) , w_awgn_matlab(1,:));
end
%%%%%%%%%%%%
clf
hold on
xlim([-4, 4]);
ylim([-4, 4]);
%received signal
scatter(X_(1, idx_tx),X_(2, idx_tx),3,"black",'o','DisplayName','Tx Constellation');
if cmplx
scatter(r_awgn_(1,:),r_awgn_(2,:),1,"red",'x','DisplayName','Matlab Mapping');
else
scatter(r_awgn_,zeros(size(r_awgn_)),1,"red",'x','DisplayName','Matlab Mapping');
end
drawnow
MI_awgn_(i_SNR)=air(X_, r_awgn_, idx_tx);
% The following also works but is slower
% MI_awgn(i_SNR)=mi_cg(s, r_awgn);
end
figure(4);
hold on
if isempty(I_shannon)
I_shannon=log2(1+db2pow(snr_meas_));
plot(snr_meas_, I_shannon, '-', 'DisplayName', 'log_2 (1+SNR)','Color','black');
end
if usegarcia
plot(SNR_dB_values, MI_awgn, '--', 'DisplayName', ['GARCIA ',num2str(M_PAM) '-QAM, AWGN']);
end
cols = linspecer(6);
if complex_constellation
plot(snr_meas_, MI_awgn_, ':', 'DisplayName', ['',num2str(M_PAM) '-QAM, AWGN'],'LineWidth',1,'Color',cols(cnt,:));
else
plot(snr_meas_, MI_awgn_, '-', 'DisplayName', ['',num2str(M_PAM) '-PAM, AWGN'],'LineWidth',1,'Color',cols(cnt,:));
end
hold off;
ylim([0,8]);
xlim([min(SNR_dB_values) max(SNR_dB_values)])
xlabel('SNR (dB)'); ylabel('Achievable rate (bits/complex dimension)');
title(['Achievable rate of ' num2str(M_PAM) '-QAM in AWGN']);
legend('Location', 'NorthWest');
yline(log2(M_PAM),'LineStyle',':','HandleVisibility','off','Color','black');
autoArrangeFigures;
cnt = cnt+1;
end
end

10
projects/ECOC_2025/theory/analytic_mpi_evaluation.m
View File

@@ -61,8 +61,10 @@ coherence_length_multiples = 0.5:0.5:ceil(L(end)/L_c);
figure();
hold on;
plot(L, avg_of_mc_variances, 'LineWidth',2, 'DisplayName','Simulation','Color',cols(1,:),'LineStyle','-');
errorbar(L, avg_of_mc_variances,std_of_mc_variances, 'LineWidth',0.7,'LineStyle','none', 'DisplayName','Simulation','Color',cols(1,:),'HandleVisibility','off');
[hl, hp] = boundedline(L, avg_of_mc_variances, std_of_mc_variances, 'alpha', 'cmap', cols(1,:));
set(hl, 'LineWidth', 2, 'DisplayName', 'Simulation');
set(hp, 'HandleVisibility', 'off', 'FaceAlpha', 0.8); % Hide patch from legend to match original behavior
plot(L_, analytic_variance, 'LineWidth',2, 'DisplayName','Analytic','Color',cols(2,:),'LineStyle','-');
xticks(coherence_length_multiples.*L_c);
@@ -76,10 +78,12 @@ else
xlabel('Interference Delay [m]', 'FontSize',12);
end
xline(L_c.*coherence_length_multiples, 'LineWidth',1.5, 'DisplayName','Coh. Length','HandleVisibility','off','Color',[0.7,0.7,0.7],'LineStyle','-');
%xline(L_c.*coherence_length_multiples, 'LineWidth',1.5,'HandleVisibility','off','Color',[0.7,0.7,0.7],'LineStyle','-');
xlim([0,L(end)]);
yline(var_sat, '-.k','LineWidth',1.5, 'DisplayName','Saturation: 2$\alpha ^2$');
grid on;
ylabel('Intensity Variance', 'FontSize',12);
title(sprintf('MPI Variance; %d MHz; SIR: %d dB',df.*1e-6,SIR_dB), 'FontSize',14);
legend('Location','southeast');
% mat2tikz_improved("C:\Users\Silas\Documents\6971e0b65b380ca6d71c837f\02_IMDD_System\tikz\mpi\analytical_mpi_variance2.tikz");

1
projects/ECOC_2025/theory/coherence_length_plot.m
View File

@@ -18,6 +18,7 @@ xlabel('Laser linewidth [MHz]','FontSize',12,'Interpreter','latex');
ylabel('Coherence length [m]','FontSize',12,'Interpreter','latex');
title('Coherence Length vs. Laser Linewidth','FontSize',14,'Interpreter','latex');
mat2tikz_improved("C:\Users\Silas\Documents\6971e0b65b380ca6d71c837f\02_IMDD_System\tikz\mpi\laser_linewidth_vs_coherence.tikz");
%% Annotate some key points
% hold on;
% freqs = [150e3, 1e6, 10e6, 50e6]; % [Hz]