diff --git a/Classes/00_signals/Signal.m b/Classes/00_signals/Signal.m index f56c09d..0a45e32 100644 --- a/Classes/00_signals/Signal.m +++ b/Classes/00_signals/Signal.m @@ -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]) diff --git a/Functions/Theory/Dissertation/dispersion_around_zdw.m b/Functions/Theory/Dissertation/dispersion_around_zdw.m index 4924f01..2c3bbf0 100644 --- a/Functions/Theory/Dissertation/dispersion_around_zdw.m +++ b/Functions/Theory/Dissertation/dispersion_around_zdw.m @@ -153,34 +153,34 @@ 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; + % 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 diff --git a/Functions/Theory/Dissertation/pmd_vs_length.m b/Functions/Theory/Dissertation/pmd_vs_length.m new file mode 100644 index 0000000..ef43b3b --- /dev/null +++ b/Functions/Theory/Dissertation/pmd_vs_length.m @@ -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') diff --git a/Libs/mutual information rate TUM/LICENSE.txt b/Libs/mutual information rate TUM/LICENSE.txt new file mode 100644 index 0000000..392b8a7 --- /dev/null +++ b/Libs/mutual information rate TUM/LICENSE.txt @@ -0,0 +1,25 @@ +Copyright (c) 2019 Francisco Javier Garcia-Gomez +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. diff --git a/Libs/mutual information rate TUM/README.md b/Libs/mutual information rate TUM/README.md new file mode 100644 index 0000000..0392605 --- /dev/null +++ b/Libs/mutual information rate TUM/README.md @@ -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. diff --git a/Libs/mutual information rate TUM/air.m b/Libs/mutual information rate TUM/air.m new file mode 100644 index 0000000..951192f --- /dev/null +++ b/Libs/mutual information rate TUM/air.m @@ -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 + diff --git a/Libs/mutual information rate TUM/example_mi_cg_complex_16qam.m b/Libs/mutual information rate TUM/example_mi_cg_complex_16qam.m new file mode 100644 index 0000000..51185e2 --- /dev/null +++ b/Libs/mutual information rate TUM/example_mi_cg_complex_16qam.m @@ -0,0 +1,124 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Copyright (c) 2019 Francisco Javier Garcia-Gomez +% 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 +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +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'); diff --git a/Libs/mutual information rate TUM/example_mi_cg_real_2D_awgn_vs_phasenoise.m b/Libs/mutual information rate TUM/example_mi_cg_real_2D_awgn_vs_phasenoise.m new file mode 100644 index 0000000..c473bdf --- /dev/null +++ b/Libs/mutual information rate TUM/example_mi_cg_real_2D_awgn_vs_phasenoise.m @@ -0,0 +1,140 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Copyright (c) 2019 Francisco Javier Garcia-Gomez +% 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 +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +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'); diff --git a/Libs/mutual information rate TUM/example_silas.m b/Libs/mutual information rate TUM/example_silas.m new file mode 100644 index 0000000..ff2faa7 --- /dev/null +++ b/Libs/mutual information rate TUM/example_silas.m @@ -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 \ No newline at end of file diff --git a/Libs/mutual information rate TUM/mi_cg.m b/Libs/mutual information rate TUM/mi_cg.m new file mode 100644 index 0000000..f0cf9f3 --- /dev/null +++ b/Libs/mutual information rate TUM/mi_cg.m @@ -0,0 +1,104 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Copyright (c) 2019 Francisco Javier Garcia-Gomez +% 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 +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +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'); diff --git a/Libs/mutual information rate TUM/silas_example_mi_cg_pam.m b/Libs/mutual information rate TUM/silas_example_mi_cg_pam.m new file mode 100644 index 0000000..e15d857 --- /dev/null +++ b/Libs/mutual information rate TUM/silas_example_mi_cg_pam.m @@ -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 \ No newline at end of file diff --git a/projects/ECOC_2025/theory/analytic_mpi_evaluation.m b/projects/ECOC_2025/theory/analytic_mpi_evaluation.m index 460f24f..56b0391 100644 --- a/projects/ECOC_2025/theory/analytic_mpi_evaluation.m +++ b/projects/ECOC_2025/theory/analytic_mpi_evaluation.m @@ -7,7 +7,7 @@ alpha = 10^(-SIR_dB/20); % Interference attenuation [linear] n_fiber = 1.467; % Refractive index c = physconst('lightspeed'); % [m/s] -L = linspace(0,250,50); % Interference delay [m] +L = linspace(0,250,50); % Interference delay [m] tau = n_fiber./c.*L; % Interference time (= tau) [s] tau_c = 1/(pi*df); % laser coherence time [s] @@ -25,21 +25,21 @@ phase_noise_std = sqrt(2*pi*df/fs); % standard dev. phase noise num_realizations = 50; % number of parallel runs monte_carlo_variance = zeros(num_realizations, length(L)); parfor r = 1:num_realizations - + % generate a realization of phase noise random walk dphi = phase_noise_std * randn(1, N + max_delay_samples); % matlab randn process has std = 1 phi = cumsum(dphi); phi_direct = phi(max_delay_samples+1 : max_delay_samples+N); var_k = zeros(1, length(L)); for t = 1:length(tau) - + nd = round( tau(t)*fs ); % delay in samples for current interference time phi_delayed = phi(max_delay_samples+1-nd : max_delay_samples+N-nd); %cut out interfering signal part (was earlier) - + E = exp(1j*phi_direct) + alpha*exp(1j*phi_delayed); % E-fields combined I = abs(E).^2; % photo current as magnitude square of E-field var_k(t) = var(I); - + end monte_carlo_variance(r, :) = var_k; end @@ -48,21 +48,23 @@ avg_of_mc_variances = mean(monte_carlo_variance, 1); std_of_mc_variances = std(monte_carlo_variance, 0, 1); %% Analytic variance -L_ = linspace(0,250,500); % Interference delay [m] +L_ = linspace(0,250,500); % Interference delay [m] tau_ = n_fiber./c.*L_; analytic_variance = 2*alpha^2 * (1 - exp(-2*pi*df.*tau_)).^2; -%% Plot +%% Plot cols = [0.3467 0.5360 0.6907 - 0.9153 0.2816 0.2878 - 0.4416 0.7490 0.4322]; + 0.9153 0.2816 0.2878 + 0.4416 0.7490 0.4322]; 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"); \ No newline at end of file diff --git a/projects/ECOC_2025/theory/coherence_length_plot.m b/projects/ECOC_2025/theory/coherence_length_plot.m index b793687..9d8fcc8 100644 --- a/projects/ECOC_2025/theory/coherence_length_plot.m +++ b/projects/ECOC_2025/theory/coherence_length_plot.m @@ -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] diff --git a/test/matched_filter_minimal.m b/test/matched_filter_minimal.m index 3f2ff9f..a7a2a84 100644 --- a/test/matched_filter_minimal.m +++ b/test/matched_filter_minimal.m @@ -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); -Digi_sig = Pform.process(symbols); -Digi_sig.spectrum("displayname",'Signal after pluse shaping','fignum',1); + +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.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);