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

Ordnerstruktur

Browse Source
This commit is contained in:
Silas Oettinghaus
2026-04-23 16:16:58 +02:00
parent 8fd68c0a54
commit 456fd02186
13 changed files with 1117 additions and 194 deletions
Download Patch File Download Diff File Expand all files Collapse all files

154
Classes/00_signals/Signal.m
View File

@@ -344,14 +344,16 @@ classdef Signal
obj
options.fignum = 2025
options.displayname = "";
options.color = [];
options.linestyle = '-';
options.normalizeToNyquist = 0;
options.addDCoffset = 0;
options.normalizeToDC = 0;
options.normalizeTo0dB = 0;
options.max_num_lines = []; % Leave empty or omit to disable line rotation
options.fft_length = [];
options.color = [];
options.linestyle = '-';
options.normalizeToNyquist = 0;
options.normalizeToSamplingRate = 0;
options.addDCoffset = 0;
options.normalizeToDC = 0;
options.normalizeTo0dB = 0;
options.show_onesided = false;
options.max_num_lines = []; % Leave empty or omit to disable line rotation
options.fft_length = [];
% --- NEW options ---
options.useWavelengthAxis (1,1) logical = false % plot x-axis in wavelength
options.lambda0_nm (1,1) double = 1310 % center wavelength [nm]
@@ -362,18 +364,21 @@ classdef Signal
end
if options.normalizeToNyquist == 0
[p_lin,f_Hz] = pwelch(obj.signal, hanning(options.fft_length), ...
options.fft_length/2, options.fft_length, ...
obj.fs, "centered", "power", "mean");
f_GHz = f_Hz*1e-9; % keep frequency vector for frequency axis
else
[p_lin,f_rad] = pwelch(obj.signal, hanning(options.fft_length), ...
options.fft_length/2, options.fft_length, ...
"centered", "power", "mean");
% In normalized mode, pwelch returns rad/sample centered on 0.
% We'll keep f_rad for the x-axis in that mode.
end
useSamplingRateAxis = options.normalizeToSamplingRate ~= 0;
useRadPerSampleAxis = options.normalizeToNyquist ~= 0 && ~useSamplingRateAxis;
if ~useRadPerSampleAxis && ~useSamplingRateAxis
[p_lin,f_Hz] = pwelch(obj.signal, hanning(options.fft_length), ...
options.fft_length/2, options.fft_length, ...
obj.fs, "centered", "power", "mean");
f_GHz = f_Hz*1e-9; % keep frequency vector for frequency axis
else
[p_lin,f_rad] = pwelch(obj.signal, hanning(options.fft_length), ...
options.fft_length/2, options.fft_length, ...
"centered", "power", "mean");
% In normalized modes, pwelch returns rad/sample centered on 0.
% Divide by 2*pi for the f/fs axis where Nyquist is 0.5.
end
% p_lin = movmean(p_lin,4);
@@ -387,45 +392,60 @@ classdef Signal
end
% --- If requested, build wavelength axis from frequency offset ---
if options.useWavelengthAxis && options.normalizeToNyquist == 0
c = physconst('LightSpeed'); % [m/s]
lambda0_m = options.lambda0_nm*1e-9; % center wavelength [m]
f_c = c / lambda0_m; % carrier frequency [Hz]
if options.useWavelengthAxis && ~useRadPerSampleAxis && ~useSamplingRateAxis
c = physconst('LightSpeed'); % [m/s]
lambda0_m = options.lambda0_nm*1e-9; % center wavelength [m]
f_c = c / lambda0_m; % carrier frequency [Hz]
% exact mapping
f_abs = f_c + f_Hz; % absolute frequency [Hz]
lambda_m = c ./ f_abs; % wavelength [m]
lambda_nm = lambda_m * 1e9; % wavelength [nm]
% assign axis
x_vec = lambda_nm(:);
x_label = "Wavelength [nm]";
% Sort to ensure axis is ascending
[x_vec, sortIdx] = sort(x_vec, 'ascend');
p_dbm = p_dbm(sortIdx, :);
else
% Frequency or normalized axes
if options.normalizeToNyquist == 0
x_vec = f_GHz;
x_label = "Frequency in GHz";
else
x_vec = f_rad; % normalized frequency in rad/sample
x_label = "Normalized Frequency";
end
end
figure(options.fignum);
ax = gca;
% assign axis
x_vec = lambda_nm(:);
x_label = "Wavelength [nm]";
dc_axis = f_Hz;
% Sort to ensure axis is ascending
[x_vec, sortIdx] = sort(x_vec, 'ascend');
p_dbm = p_dbm(sortIdx, :);
dc_axis = dc_axis(sortIdx);
else
% Frequency or normalized axes
if ~useRadPerSampleAxis && ~useSamplingRateAxis
x_vec = f_GHz;
x_label = "Frequency in GHz";
dc_axis = f_GHz;
elseif useSamplingRateAxis
x_vec = f_rad ./ (2*pi);
x_label = "Normalized Frequency f/fs";
dc_axis = x_vec;
else
x_vec = f_rad; % normalized frequency in rad/sample
x_label = "Normalized Frequency [rad/sample]";
dc_axis = x_vec;
end
end
if options.show_onesided
keep_idx = dc_axis >= 0;
x_vec = x_vec(keep_idx);
dc_axis = dc_axis(keep_idx);
p_dbm = p_dbm(keep_idx, :);
end
figure(options.fignum);
ax = gca;
hold on
p_dbm = p_dbm+options.addDCoffset;
if options.normalizeToDC
[~,min_idx]=min(abs(f_GHz));
pow_at_dc = p_dbm(min_idx);
p_dbm = p_dbm-pow_at_dc;
end
p_dbm = p_dbm+options.addDCoffset;
if options.normalizeToDC
[~,min_idx]=min(abs(dc_axis));
pow_at_dc = p_dbm(min_idx);
p_dbm = p_dbm-pow_at_dc;
end
% p_dbm = movmean(p_dbm,10);
for s = 1:min(size(p_dbm))
@@ -448,17 +468,27 @@ classdef Signal
% Axis labels and limits
xlabel(x_label);
if options.useWavelengthAxis && options.normalizeToNyquist == 0
xlim([min(x_vec) max(x_vec)]);
else
if options.normalizeToNyquist == 0
% Keep your existing freq handling (you can fine-tune as needed)
% xlim([-128 128]); % example for 256 GSa/s if desired
xlim([min(x_vec) max(x_vec)]);
else
xlim([-pi, pi]);
end
end
if options.useWavelengthAxis && ~useRadPerSampleAxis && ~useSamplingRateAxis
xlim([min(x_vec) max(x_vec)]);
else
if ~useRadPerSampleAxis && ~useSamplingRateAxis
% Keep your existing freq handling (you can fine-tune as needed)
% xlim([-128 128]); % example for 256 GSa/s if desired
xlim([min(x_vec) max(x_vec)]);
elseif useSamplingRateAxis
if options.show_onesided
xlim([0, 0.5]);
else
xlim([-0.5, 0.5]);
end
else
if options.show_onesided
xlim([0, pi]);
else
xlim([-pi, pi]);
end
end
end
ylabel(ylab);

64
Classes/01_transmit/Signalgenerator.m
View File

@@ -21,21 +21,22 @@ classdef Signalgenerator
% "dimension", 2, ...
% "randkey", 3).process();
%
% % PRMS bit matrix for PAM-4 mapping. If length is omitted, the
% % output length is derived from order, as in the old PAMsource.
% % PRMS bit matrix for PAM-4 mapping. M selects the mapper-ready bit
% % shape. If length is omitted, the output length is derived from
% % order, as in the old PAMsource.
% bits = Signalgenerator( ...
% "form", signalform.prms, ...
% "dimension", 2, ...
% "M", 4, ...
% "order", 7).process();
% symbols = PAMmapper(4, 0).map(bits);
%
% % PRMS for PAM-6: 5 bits map to 2 PAM-6 symbols
% % PRMS for PAM-6: 5 bits map to 2 PAM-6 symbols. M = 6 returns
% % the 1-D bit vector expected by PAMmapper.
% bits = Signalgenerator( ...
% "form", signalform.prms, ...
% "dimension", 5, ...
% "M", 6, ...
% "order", 7).process();
% bitVector = reshape(bits.signal.', [], 1);
% symbols = PAMmapper(6, 0).map(bitVector);
% symbols = PAMmapper(6, 0).map(bits);
%
% % Map, demap, and check BER
% mapper = PAMmapper(8, 0);
@@ -54,6 +55,7 @@ classdef Signalgenerator
fs
fsig
dimension
M
order
randkey
skip
@@ -72,6 +74,7 @@ classdef Signalgenerator
options.fs double = 1000 %Hz sampling
options.fsig double = 50 % Hz fundamental frex e.g. of the sine or sawtooth
options.dimension double = 1
options.M double = []
options.order double = 7
options.randkey double = 0
options.skip double = 0
@@ -86,6 +89,11 @@ classdef Signalgenerator
end
end
if ~isempty(obj.M)
obj.validate_pam_format();
obj.dimension = obj.mapper_bit_dimension();
end
if isempty(obj.length)
obj.length = obj.default_length();
end
@@ -136,7 +144,8 @@ classdef Signalgenerator
case signalform.prms
signal = obj.build_prms_data();
end
signal = obj.format_for_mapper(signal);
end
@@ -169,12 +178,15 @@ classdef Signalgenerator
% PAM-8: dimension = 3
%
% PAM-6 is special in this codebase: the mapper consumes 5 bits and
% maps them to two PAM-6 symbols. Use dimension = 5, then flatten
% the output before calling PAMmapper(6, ...).map(...).
% maps them to two PAM-6 symbols. Prefer M = 6 to generate the 1-D
% vector expected by PAMmapper(6, ...).map(...). The legacy
% dimension = 5 path still returns the unflattened bit matrix.
%
% Relevant Signalgenerator options:
% length Number of generated PRMS bit groups / rows.
% dimension Number of parallel bit streams per row.
% M Optional PAM format. When set, dimension is selected
% automatically and PAM-6 output is mapper-ready.
% order User-facing sequence order, matching the old
% PAMsource behavior. If length is omitted, it sets the
% output length. Internally it is converted to the PRMS
@@ -252,10 +264,17 @@ classdef Signalgenerator
end
end
function validate_pam_format(obj)
if ~ismember(obj.M, [2 4 6 8 16])
error("Signalgenerator:InvalidPAMFormat", ...
"M must be one of 2, 4, 6, 8, or 16.");
end
end
function length = default_length(obj)
switch obj.form
case {signalform.random, signalform.prms}
if obj.dimension == 5
if obj.is_pam6_shape()
length = 2^max(0, obj.order - 2);
else
length = 2^max(0, obj.order - 1);
@@ -266,7 +285,9 @@ classdef Signalgenerator
end
function prms_order = internal_prms_order(obj)
if obj.dimension == 5
if ~isempty(obj.M)
bits_per_symbol = log2(obj.M);
elseif obj.dimension == 5
bits_per_symbol = log2(6);
else
bits_per_symbol = obj.dimension;
@@ -275,6 +296,25 @@ classdef Signalgenerator
prms_order = max(1, floor(obj.order / bits_per_symbol));
end
function dimension = mapper_bit_dimension(obj)
if obj.M == 6
dimension = 5;
else
dimension = log2(obj.M);
end
end
function tf = is_pam6_shape(obj)
tf = (~isempty(obj.M) && obj.M == 6) || (isempty(obj.M) && obj.dimension == 5);
end
function signal = format_for_mapper(obj, signal)
if ~isempty(obj.M) && obj.M == 6
signal = reshape(signal.', [], 1);
signal = signal(1:end - mod(numel(signal), 5));
end
end
end

36
Functions/EQ_visuals/showErrorBurstCount.m
View File

@@ -4,6 +4,7 @@ arguments
ref_symbols
options.fignum (1,1) double = NaN % Default to NaN if not provided
options.displayname (1,:) char = '' % Default to an empty string if not provided
options.max_burst_length (1,1) double = 20
end
% Determine the figure number to use or create a new figure
@@ -18,30 +19,17 @@ end
eq_signal = PAMmapper(M,0).quantize(eq_signal);
end
diff_indices = eq_signal.signal == ref_symbols.signal;
errpos = find(eq_signal.signal(:) ~= ref_symbols.signal(:));
burst_len = 1:options.max_burst_length;
burst_count = count_error_bursts(errpos, options.max_burst_length);
% Identify the start of new sequences (when the difference is not 1)
sequence_starts = [1, find(diff_indices ~= 1) + 1]; % Include the first index
sequence_ends = [sequence_starts(2:end) - 1, length(errpos)]; % Calculate end indices
% Initialize burst count and print bursts longer than 10
burst_len = 1:10;
burst_count = zeros(length(burst_len),1);
for t = 1:numel(burst_len)
for i = 1:length(sequence_starts)
% Extract current sequence
current_burst = errpos(sequence_starts(i):sequence_ends(i));
% Check if the sequence length matches criterion
if length(current_burst) == burst_len(t)
burst_count(t) = burst_count(t) + 1;
end
end
bar(burst_len, burst_count, "DisplayName", options.displayname);
grid on;
xlabel("Error burst length");
ylabel("Count");
title("Error burst count");
if ~isempty(options.displayname)
legend("show");
end
burst_symbols = burst_count .* burst_len';
burst_rate = burst_symbols ;%./ length(rx_signal);
end
end

6
Functions/Metrics/count_error_bursts.m
View File

@@ -3,10 +3,14 @@ function burst_count = count_error_bursts(err_pos, max_burst_length)
% max_burst_length: Maximum length for which bursts are counted (e.g., 10)
% Sort the error positions to ensure they're in increasing order
err_pos = sort(err_pos);
err_pos = sort(err_pos(:).');
% Initialize burst_count array to hold the counts for each burst length
burst_count = zeros(1, max_burst_length);
if isempty(err_pos)
return
end
% Find the differences between consecutive error positions
diffs = diff(err_pos);

47
Tests/00_signals/Signal_test.m
View File

@@ -87,6 +87,53 @@ classdef Signal_test < IMDDTestCase
testCase.verifyEqual(height(sig.logbook), 1);
end
function spectrumCanUseSamplingRateNormalizedFrequencyAxis(testCase)
fig = figure("Visible", "off");
cleanup = onCleanup(@() close(fig));
fignum = fig.Number;
sig = Signal(randn(4096, 1), "fs", 8e9);
sig.spectrum( ...
"fignum", fignum, ...
"normalizeToSamplingRate", 1, ...
"normalizeTo0dB", 1, ...
"fft_length", 512);
ax = findobj(fig, "Type", "Axes");
lineObj = findobj(ax, "Type", "Line");
testCase.verifyEqual(ax.XLim, [-0.5 0.5], "AbsTol", 1e-12);
testCase.verifyEqual(string(ax.XLabel.String), "Normalized Frequency f/fs");
testCase.verifyGreaterThanOrEqual(min(lineObj.XData), -0.5);
testCase.verifyLessThanOrEqual(max(lineObj.XData), 0.5);
testCase.verifyLessThanOrEqual(abs(min(lineObj.XData) + 0.5), 1/512);
testCase.verifyLessThanOrEqual(abs(max(lineObj.XData) - 0.5), 1/512);
end
function spectrumCanShowOneSidedSamplingRateNormalizedAxis(testCase)
fig = figure("Visible", "off");
cleanup = onCleanup(@() close(fig));
fignum = fig.Number;
sig = Signal(randn(4096, 1), "fs", 8e9);
sig.spectrum( ...
"fignum", fignum, ...
"normalizeToSamplingRate", 1, ...
"show_onesided", true, ...
"normalizeTo0dB", 1, ...
"fft_length", 512);
ax = findobj(fig, "Type", "Axes");
lineObj = findobj(ax, "Type", "Line");
testCase.verifyEqual(ax.XLim, [0 0.5], "AbsTol", 1e-12);
testCase.verifyGreaterThanOrEqual(min(lineObj.XData), 0);
testCase.verifyLessThanOrEqual(max(lineObj.XData), 0.5);
testCase.verifyLessThanOrEqual(abs(max(lineObj.XData) - 0.5), 1/512);
end
function castingRoundTripPreservesSignalThroughAllDomains(testCase)
infoIn = Informationsignal([0; 1; 1; 0]);

29
Tests/01_transmit/Signalgenerator_test.m
View File

@@ -8,6 +8,7 @@ classdef Signalgenerator_test < IMDDTestCase
testCase.verifyEqual(gen.fs, 1000);
testCase.verifyEqual(gen.fsig, 50);
testCase.verifyEqual(gen.dimension, 1);
testCase.verifyEmpty(gen.M);
testCase.verifyEqual(gen.order, 7);
testCase.verifyEqual(gen.randkey, 0);
testCase.verifyEqual(gen.skip, 0);
@@ -102,6 +103,34 @@ classdef Signalgenerator_test < IMDDTestCase
testCase.verifyEqual(pam6Gen.length, 2^(7 - 2));
end
function pamFormatSelectsMapperReadyDimension(testCase)
pam4Gen = Signalgenerator( ...
"form", signalform.prms, ...
"M", 4, ...
"order", 7);
pam6Gen = Signalgenerator( ...
"form", signalform.prms, ...
"M", 6, ...
"order", 7);
testCase.verifyEqual(pam4Gen.dimension, 2);
testCase.verifyEqual(pam4Gen.length, 2^(7 - 1));
testCase.verifyEqual(pam6Gen.dimension, 5);
testCase.verifyEqual(pam6Gen.length, 2^(7 - 2));
end
function pam6FormatReturnsMapperReadyVector(testCase)
gen = Signalgenerator( ...
"form", signalform.prms, ...
"M", 6, ...
"order", 7);
signal = gen.build_signal();
testCase.verifySize(signal, [5 * 2^(7 - 2), 1]);
testCase.verifyEqual(mod(numel(signal), 5), 0);
end
function processWrapsPrmsSignalInInformationsignal(testCase)
gen = Signalgenerator( ...
"form", signalform.prms, ...

18
Theory/Mapping_Coding/partial_response/duobinary_back2back.m Normal file
View File

@@ -0,0 +1,18 @@
for M = [2 4 8]
bits = Signalgenerator("form", signalform.prms,"M", M,"order", 15).process();
symbols = PAMmapper(M,0).map(bits);
symbols_pre = Duobinary().precode(symbols);
symbols_db = Duobinary().encode(symbols_pre);
symbols_rx = Duobinary().decode(symbols_db);
% Vergleichen von b(n) und d_dec(n)
bits_rx = PAMmapper(M,0).demap(symbols_rx);
[~,~,ber,~] = calc_ber(bits.signal,bits_rx.signal,"skip_front",10,"skip_end",10,"returnErrorLocation",1);
disp(['BER: ',sprintf('%.1E',ber),' - - PAM-',num2str(M)]);
end

26
Theory/Mapping_Coding/partial_response/duobinary_minimal_example.m
View File

@@ -9,8 +9,7 @@ Tx_bits = Signalgenerator( ...
"M", M, ...
"order", 15).process();
%%%%% Duobinary %%%%%%
close all
Symbols_tx = PAMmapper(M,0,"eth_style",0).map(Tx_bits);
Symbols_tx.fs = fsym;
@@ -29,8 +28,12 @@ end
for n = 10
% Start from the transmitted symbol sequence for this impairment case.
Symbols_rx = Symbols_tx;
% Insert n deterministic symbol errors near the beginning of the sequence.
% Each replacement is chosen from another symbol position with a different
% amplitude, so this produces actual symbol decisions errors.
pos = 1;
if n~=0
for pos = 1:n
@@ -45,9 +48,17 @@ for n = 10
end
end
% Keep track of the positions where the artificial symbol errors were
% inserted before duobinary encoding/decoding is applied.
error_positions = ~(Symbols_rx.signal == Symbols_tx.signal);
error_positions = find(error_positions==1);
% Receiver-side duobinary processing depends on what was already done at
% the transmitter side:
% - db_precoded: the sequence is precoded only, so emulate the channel PR
% encoder first and then apply memoryless duobinary decoding.
% - db_encoded: the sequence already contains duobinary levels, so only
% the memoryless duobinary decoder is applied.
switch precode
case db_mode.db_precoded
@@ -58,9 +69,12 @@ for n = 10
end
% Convert the recovered PAM symbols back to bits with the framework mapper.
Rx_bits = PAMmapper(M,0).demap(Symbols_rx);
%%%%% Check BER of Bit Sequence %%%%%
% Skip a few edge samples because the first/last symbols can include
% state/transient effects from the duobinary recursion.
[~,error_num(n+1),ber,error_pos] = calc_ber(Tx_bits.signal,Rx_bits.signal,"skip_front",10,"skip_end",10,"returnErrorLocation",1);
% disp(['BER: ',sprintf('%.1E',ber),sprintf(' - Num. Err: %.1d',error_num(n+1)-2),' - - PAM-',num2str(M)]);
@@ -70,9 +84,12 @@ end
% Plot bit-level and symbol-level comparisons for a quick visual sanity check.
figure(3);
clf
%sgtitle(['BER: ',num2str(ber),' // Error is at position: ',num2str(error_pos),''])
% Compare the first few transmitted and received bits.
subplot(2,2,1)
hold on
title('First Bits')
@@ -80,6 +97,7 @@ stairs(Tx_bits.signal(100:150,1),'LineStyle','-','LineWidth',2,'DisplayName','Tx
stairs(Rx_bits.signal(100:150,1),'LineWidth',2,'DisplayName','Rx Bits','LineStyle',':')
legend
% Compare the last few transmitted and received bits.
subplot(2,2,2)
title('Last Bits')
hold on
@@ -87,6 +105,7 @@ stairs(Tx_bits.signal(end-50:end,1),'LineStyle','-','LineWidth',2,'DisplayName',
stairs(Rx_bits.signal(end-50:end,1),'LineWidth',2,'DisplayName','Rx Bits','LineStyle',':')
legend
% Compare the first few symbol decisions after the duobinary processing chain.
subplot(2,2,3)
hold on
title('First Symbols Compare')
@@ -94,9 +113,10 @@ stairs(Symbols_tx.signal(1:100,1),'LineWidth',2,'DisplayName','Tx Symbols','Line
stairs(Symbols_rx.signal(1:100,1),'LineStyle',':','LineWidth',2,'DisplayName','Rx Symbols');
legend
% Compare the last few symbol decisions after the duobinary processing chain.
subplot(2,2,4)
hold on
title('Last Symbols Compare')
stairs(Symbols_tx.signal(end-50:end,1),'LineWidth',2,'DisplayName','Tx Symbols','LineStyle','-')
stairs(Symbols_rx.signal(end-50:end,1),'LineStyle',':','LineWidth',2,'DisplayName','Rx Symbols');
legend
legend

73
Theory/Mapping_Coding/partial_response/pam_6_differential_code_understand.m
View File

@@ -1,73 +0,0 @@
M = 6;
data = [1,2,3,4,5,6];
M = 6;
bitpattern = [];
s = RandStream('twister','Seed',1);
for i = 1:log2(M)
N = 2^(12-1); %length of prbs
bitpattern(:,i) = randi(s,[0 1], N, 1);
end
if M == 6
bitpattern = reshape(bitpattern',[],1);
bitpattern = bitpattern(1:end-mod(length(bitpattern),5));
end
bits = Informationsignal(bitpattern);
symbols = PAMmapper(M,0).map(bits);
symbols_tx_prec = Duobinary().precode(symbols);
% all possible transitions (for now 36, including the "edges"
% of the QAM 32 constellation)
states = PAMmapper(6,0,"eth_style",0).levels;
pam6transitions = combvec(states,states)'; % pam6transitions =
% [-5 -5;
% -3 -5;
% -1 -5; ...
pam6transitions_serial = reshape(pam6transitions',[],1);
data = pam6transitions_serial;
data = round(data);
b = min(data);
data = data - b;
data = data ./ 2;
% THIS WAS USED!
bk = zeros(size(data));
for k = 2:numel(data)
bk(k) = mod(data(k)-bk(k-1),M);
end
%% State Analysis
x = bk;%symbols_tx_prec.signal;
levels = sort(unique(x)).'; % or provide known 1x6 level values
[~,ix] = min(abs(x - levels),[],2);
x = levels(ix); % snapped/quantized
%% TRANSITION COUNTS & PROBABILITIES
K = numel(levels);
% map to state indices 1..K
[tf, idx] = ismember(x, levels);
idx = idx(:);
from = idx(1:end-1);
to = idx(2:end);
from = idx(1:2:end);
to = idx(2:2:end);
% counts C(from,to)
C = accumarray([from,to], 1, [K K], @sum, 0);
% row-stochastic transition matrix P(to|from)
rowSums = sum(C,2);
P = C ./ max(rowSums,1);
%% 1) HEATMAP (which transitions are more probable?)
figure('Name','Transition Probabilities (to | from)');
h = heatmap(levels, levels, P, 'Colormap', parula, 'ColorbarVisible','on');
colormap(gca,[[1,1,1];flip(cbrewer2('Spectral',100))]);clim([0,ceil(max(P(:))*10)/10]);
h.XLabel = 'From state (level)';
h.YLabel = 'To state (level)';
h.Title = 'P(to | from)';

835
Theory/Mapping_Coding/partial_response/partial_response_precoder.m Normal file
View File

@@ -0,0 +1,835 @@
% Generalized duobinary partial-response example for PAM-M.
% Symbol-domain implementation for F(D) = (1 + D)^n, M in {2, 4, 8}.
clear; clc;
repoRoot = fileparts(fileparts(fileparts(mfilename("fullpath"))));
addpath(genpath(repoRoot));
rng(7);
M_list = [2,4,8];
pr_orders = [1,2,3];
Nsym = 100000;
userOrder = 17;
emulate_precoder = 0;
burst_lengths = 1:50;
injected_error_count = 100;
plot_burst_analysis = 1;
burst_randkey = 140;
burst_min_gap = 100;
fprintf('Generalized duobinary PR burst-resilience demo\n');
fprintf('emulate_precoder = %d\n', emulate_precoder);
fprintf('Injected precoder-symbol errors per burst length = %d\n', injected_error_count);
fprintf('Burst randkey = %d, minimum gap = %d symbols\n', burst_randkey, burst_min_gap);
burst_results = struct();
burst_pattern = buildCombinedBurstPattern( ...
Nsym, burst_lengths, injected_error_count, ...
"randkey", burst_randkey, ...
"min_gap", burst_min_gap);
for pr_order = pr_orders
fprintf('\nF(D) = (1 + D)^%d\n', pr_order);
for M = M_list
[txBits, txSymbols, mapper, d_tx] = generatePamSymbols(M, Nsym, userOrder);
if emulate_precoder
[~, y_amp, pr_alphabet, y_idx, a_ref] = emulatePrecoder(d_tx, M, [], pr_order);
refBits = demapSymbolIndices(a_ref, txSymbols, mapper, M);
else
a_ref = d_tx;
refBits = txBits;
u = precode(a_ref, M, [], pr_order);
[y_amp, pr_alphabet, y_idx] = encode(u, M, [], pr_order);
end
a_hat0 = memorylessdecode(y_amp, M, [], pr_order);
rxBits0 = demapSymbolIndices(a_hat0, txSymbols, mapper, M);
[~, errors0, ber0] = calc_ber(rxBits0.signal, refBits.signal);
assert(isequal(a_hat0, a_ref), ...
'Noiseless reconstruction failed for M=%d, order=%d.', M, pr_order);
assert(errors0 == 0 && ber0 == 0, ...
'Noiseless BER check failed for M=%d, order=%d.', M, pr_order);
assert(numel(pr_alphabet) == (M-1)*2^pr_order + 1, ...
'Wrong PR alphabet size for M=%d, order=%d.', M, pr_order);
results = simulatePrecoderErrorBursts( ...
a_ref, M, pr_order, txSymbols, mapper, refBits, ...
burst_pattern);
fprintf('M=%d before encode/decode: BER = %.3g (%d/%d bits), symbol errors = %d, bursts = %d, max burst = %d\n', ...
M, results.input_ber, results.input_bit_errors, results.input_checked_bits, ...
results.input_symbol_errors, results.input_burst_total, results.input_max_burst_length);
fprintf('M=%d after encode/decode: BER = %.3g (%d/%d bits), symbol errors = %d, bursts = %d, max burst = %d\n', ...
M, results.ber, results.bit_errors, results.checked_bits, ...
results.symbol_errors, results.output_burst_total, results.output_max_burst_length);
printBurstCountReport(results, M, pr_order);
burst_results(pr_order, M).results = results;
if plot_burst_analysis
plotPrecoderBurstAnalysis(results, pr_order, M);
end
end
end
if plot_burst_analysis
plotBurstSummary(burst_results, pr_orders, M_list);
end
% runDuobinaryPartialResponseTests();
function pam_levels = getPamLevels(M)
%GETPAMLEVELS Symmetric PAM amplitudes for natural indices 0...(M-1).
validatePamOrder(M);
pam_levels = -(M-1):2:(M-1);
end
function [txBits, txSymbols, mapper, a] = generatePamSymbols(M, Nsym, userOrder)
%GENERATEPAMSYMBOLS Framework-based bit generation and PAM mapping.
validatePamOrder(M);
txBits = Signalgenerator( ...
"form", signalform.prms, ...
"M", M, ...
"length", Nsym, ...
"order", userOrder).process();
mapper = PAMmapper(M, 0);
txSymbols = mapper.map(txBits);
a = amplitudeToSymbolIndex(txSymbols.signal .* mapper.scaling, M);
end
function rxBits = demapSymbolIndices(a_hat, txSymbolsTemplate, mapper, M)
%DEMAPSYMBOLINDICES Convert decoded indices back to framework PAM symbols.
rxSymbols = txSymbolsTemplate;
rxSymbols.signal = symbolIndexToAmplitude(a_hat, M) ./ mapper.scaling;
rxBits = mapper.demap(rxSymbols);
end
function results = simulatePrecoderErrorBursts(a_ref, M, pr_order, txSymbols, mapper, refBits, burst_pattern)
%SIMULATEPRECODERERRORBURSTS Inject deterministic symbol errors before PR encoding.
u_ref = precode(a_ref, M, [], pr_order);
results = struct();
results.burst_lengths = burst_pattern.burst_lengths(:);
results.injected_error_count_per_length = burst_pattern.error_count_per_length;
results.injected_total_errors = numel(burst_pattern.err_pos);
results.injected_burst_count = count_error_bursts( ...
burst_pattern.err_pos, max(burst_pattern.burst_lengths));
[u_err, injected_err_pos] = injectSymbolErrorBursts(u_ref, M, burst_pattern);
u_ref_bits = demapSymbolIndices(u_ref, txSymbols, mapper, M);
inputBits = demapSymbolIndices(u_err, txSymbols, mapper, M);
[inputCheckedBits, inputBitErrors, inputBer] = calc_ber( ...
inputBits.signal, u_ref_bits.signal, "returnErrorLocation", 1);
input_err_pos = find(u_err(:) ~= u_ref(:));
input_burst_lengths = getErrorBurstLengths(input_err_pos);
[y_amp_err, ~, y_idx_err] = encode(u_err, M, [], pr_order);
a_hat = memorylessdecode(y_amp_err, M, [], pr_order);
rxBits = demapSymbolIndices(a_hat, txSymbols, mapper, M);
[checkedBits, bitErrors, ber,errorIndice] = calc_ber( ...
rxBits.signal, refBits.signal, "returnErrorLocation", 1);
output_err_pos = find(a_hat(:) ~= a_ref(:));
output_burst_lengths = getErrorBurstLengths(output_err_pos);
results.input_symbol_errors = numel(input_err_pos);
results.input_symbol_ser = results.input_symbol_errors / numel(a_ref);
results.input_bit_errors = inputBitErrors;
results.input_checked_bits = inputCheckedBits;
results.input_ber = inputBer;
results.input_burst_count = count_error_bursts( ...
input_err_pos, max(burst_pattern.burst_lengths));
results.input_burst_total = numel(input_burst_lengths);
results.input_max_burst_length = 0;
if ~isempty(input_burst_lengths)
results.input_max_burst_length = max(input_burst_lengths);
end
results.symbol_errors = numel(output_err_pos);
results.symbol_ser = results.symbol_errors / numel(a_ref);
results.bit_errors = bitErrors;
results.checked_bits = checkedBits;
results.ber = ber;
results.output_burst_count = count_error_bursts( ...
output_err_pos, max(burst_pattern.burst_lengths));
results.output_burst_total = numel(output_burst_lengths);
results.output_max_burst_length = 0;
if ~isempty(output_burst_lengths)
results.output_max_burst_length = max(output_burst_lengths);
end
results.injected_err_pos = injected_err_pos;
results.output_err_pos = output_err_pos;
results.last_a_hat = a_hat;
results.last_a_ref = a_ref;
results.last_y_idx_err = y_idx_err;
end
function [u_err, err_pos] = injectSymbolErrorBursts(u_ref, M, burst_pattern)
%INJECTSYMBOLERRORBURSTS Apply one precomputed burst pattern in one operation.
validateSymbolVector(u_ref, M, 'u_ref');
original_size = size(u_ref);
u_work = double(u_ref(:));
N = numel(u_work);
err_pos = burst_pattern.err_pos(:);
if any(err_pos < 1) || any(err_pos > N)
error('injectSymbolErrorBursts:InvalidPattern', ...
'burst_pattern.err_pos contains positions outside the sequence.');
end
u_err = u_work;
u_err(err_pos) = mod(u_err(err_pos) + 1, M);
u_err = reshape(u_err, original_size);
end
function burst_pattern = buildCombinedBurstPattern(N, burst_lengths, error_count_per_length, options)
%BUILDCOMBINEDBURSTPATTERN One common sequence containing all burst lengths.
arguments
N (1,1) double
burst_lengths
error_count_per_length (1,1) double
options.randkey = []
options.min_gap (1,1) double = 20
options.max_placement_attempts (1,1) double = 20000
end
if any(burst_lengths < 1) || any(burst_lengths ~= floor(burst_lengths))
error('buildCombinedBurstPattern:InvalidBurstLengths', ...
'burst_lengths must contain positive integers.');
end
if ~(isscalar(error_count_per_length) && ...
error_count_per_length == floor(error_count_per_length) && ...
error_count_per_length >= 1)
error('buildCombinedBurstPattern:InvalidErrorCount', ...
'error_count_per_length must be a positive integer.');
end
if options.min_gap < 0 || options.min_gap ~= floor(options.min_gap)
error('buildCombinedBurstPattern:InvalidGap', ...
'min_gap must be a nonnegative integer.');
end
burst_lengths = burst_lengths(:).';
burst_plan = [];
burst_label_plan = [];
for burst_length = burst_lengths
remaining = error_count_per_length;
while remaining > 0
current_len = min(burst_length, remaining);
burst_plan(end + 1, 1) = current_len; %#ok<AGROW>
burst_label_plan(end + 1, 1) = burst_length; %#ok<AGROW>
remaining = remaining - current_len;
end
end
required_span = sum(burst_plan) + max(0, numel(burst_plan) - 1) * options.min_gap;
if required_span > N
error('buildCombinedBurstPattern:SequenceTooShort', ...
'Sequence too short for the requested random burst pattern and gap.');
end
if isempty(options.randkey)
rs = [];
else
rs = RandStream('mt19937ar', 'Seed', options.randkey);
end
placement_order = randomPermutation(numel(burst_plan), rs);
burst_plan = burst_plan(placement_order);
burst_label_plan = burst_label_plan(placement_order);
burst_starts = zeros(numel(burst_plan), 1);
burst_ends = zeros(numel(burst_plan), 1);
for b = 1:numel(burst_plan)
burst_len = burst_plan(b);
placed = false;
for attempt = 1:options.max_placement_attempts
start_pos = randomInteger(1, N - burst_len + 1, rs);
end_pos = start_pos + burst_len - 1;
if isBurstPlacementFree(start_pos, end_pos, burst_starts(1:b-1), ...
burst_ends(1:b-1), options.min_gap)
burst_starts(b) = start_pos;
burst_ends(b) = end_pos;
placed = true;
break
end
end
if ~placed
error('buildCombinedBurstPattern:PlacementFailed', ...
'Could not place all bursts. Reduce min_gap or error count.');
end
end
[burst_starts, sort_idx] = sort(burst_starts);
burst_ends = burst_ends(sort_idx);
burst_label_plan = burst_label_plan(sort_idx);
err_pos = [];
err_label = [];
for b = 1:numel(burst_starts)
current_pos = (burst_starts(b):burst_ends(b)).';
err_pos = [err_pos; current_pos]; %#ok<AGROW>
err_label = [err_label; repmat(burst_label_plan(b), numel(current_pos), 1)]; %#ok<AGROW>
end
burst_pattern = struct();
burst_pattern.err_pos = err_pos;
burst_pattern.err_label = err_label;
burst_pattern.burst_lengths = burst_lengths;
burst_pattern.error_count_per_length = error_count_per_length;
burst_pattern.randkey = options.randkey;
burst_pattern.min_gap = options.min_gap;
burst_pattern.burst_starts = burst_starts;
burst_pattern.burst_ends = burst_ends;
end
function tf = isBurstPlacementFree(start_pos, end_pos, burst_starts, burst_ends, min_gap)
%ISBURSTPLACEMENTFREE True when the candidate interval is gap-separated.
if isempty(burst_starts)
tf = true;
return
end
tf = all(end_pos + min_gap < burst_starts | start_pos - min_gap > burst_ends);
end
function p = randomPermutation(n, rs)
%RANDOMPERMUTATION randperm with optional RandStream.
if isempty(rs)
p = randperm(n);
else
p = randperm(rs, n);
end
end
function x = randomInteger(lo, hi, rs)
%RANDOMINTEGER randi with optional RandStream.
if isempty(rs)
x = randi([lo hi]);
else
x = randi(rs, [lo hi]);
end
end
function burst_lengths = getErrorBurstLengths(err_pos)
%GETERRORBURSTLENGTHS Return lengths of all contiguous error bursts.
err_pos = sort(err_pos(:).');
if isempty(err_pos)
burst_lengths = [];
return
end
diffs = diff(err_pos);
burst_starts = [1, find(diffs > 1) + 1];
burst_ends = [find(diffs > 1), numel(err_pos)];
burst_lengths = burst_ends - burst_starts + 1;
end
function u = precode(a, M, init_state, pr_order)
%PRECODE Symbol-domain modulo precoder for F(D)=(1+D)^n.
% h = coefficients of (1 + D)^n
% u(k) = mod(a(k) - sum_{i=1}^n h(i+1) u(k-i), M)
if nargin < 4 || isempty(pr_order)
pr_order = 1;
end
if nargin < 3
init_state = [];
end
validateSymbolVector(a, M, 'a');
pr_order = validatePrOrder(pr_order);
state = validateInitState(init_state, M, pr_order);
original_size = size(a);
a_work = double(a(:).');
h = getPartialResponseTaps(pr_order);
u = zeros(size(a_work));
for k = 1:numel(a_work)
% Move the PR recursion to the transmitter.
feedback = sum(h(2:end) .* state);
u(k) = mod(a_work(k) - feedback, M);
state = [u(k) state(1:end-1)];
end
u = reshape(u, original_size);
end
function [y_amp, pr_alphabet, y_idx] = encode(u, M, init_state, pr_order)
%ENCODE Apply F(D)=(1+D)^n to an already precoded sequence.
if nargin < 4 || isempty(pr_order)
pr_order = 1;
end
if nargin < 3
init_state = [];
end
validateSymbolVector(u, M, 'u');
pr_order = validatePrOrder(pr_order);
state = validateInitState(init_state, M, pr_order);
original_size = size(u);
u_work = double(u(:).');
h = getPartialResponseTaps(pr_order);
pam_levels = getPamLevels(M);
pr_alphabet = getDuobinaryAlphabet(M, pr_order);
y_idx = zeros(size(u_work));
y_amp = zeros(size(u_work));
for k = 1:numel(u_work)
pr_state = [u_work(k) state];
y_idx(k) = sum(h .* pr_state);
y_amp(k) = sum(h .* pam_levels(pr_state + 1));
state = pr_state(1:end-1);
end
y_idx = reshape(y_idx, original_size);
y_amp = reshape(y_amp, original_size);
end
function [a_hat, y_hat_amp, y_hat_idx] = memorylessdecode(r, M, init_state, pr_order)
%MEMORYLESSDECODE Slice to PR alphabet, then decode modulo M.
% a_hat(k) = mod(y_hat_idx(k), M)
if nargin < 4 || isempty(pr_order)
pr_order = 1;
end
if nargin < 3
init_state = [];
end
validatePamOrder(M);
pr_order = validatePrOrder(pr_order);
validateRealVector(r, 'r');
validateInitState(init_state, M, pr_order);
original_size = size(r);
r_work = double(r(:).');
[pr_alphabet, idx_for_alphabet] = getDuobinaryAlphabet(M, pr_order);
y_hat_amp = zeros(size(r_work));
y_hat_idx = zeros(size(r_work));
for k = 1:numel(r_work)
[~, level_pos] = min(abs(r_work(k) - pr_alphabet));
y_hat_amp(k) = pr_alphabet(level_pos);
y_hat_idx(k) = idx_for_alphabet(level_pos);
end
a_hat = mod(y_hat_idx, M);
a_hat = reshape(a_hat, original_size);
y_hat_amp = reshape(y_hat_amp, original_size);
y_hat_idx = reshape(y_hat_idx, original_size);
end
function a_hat = slicePamAmplitude(r, M)
%SLICEPAMAMPLITUDE Nearest-neighbor PAM slicing to natural symbol indices.
validatePamOrder(M);
validateRealVector(r, 'r');
original_size = size(r);
r_work = double(r(:).');
pam_levels = getPamLevels(M);
a_hat = zeros(size(r_work));
for k = 1:numel(r_work)
[~, level_pos] = min(abs(r_work(k) - pam_levels));
a_hat(k) = level_pos - 1;
end
a_hat = reshape(a_hat, original_size);
end
function plotPrecoderBurstAnalysis(results, pr_order, M)
%PLOTPRECODERBURSTANALYSIS Compact BER and burst-count visualization.
figure(300 + 10*pr_order + M);
clf;
tiledlayout(2, 1);
nexttile;
bar(results.burst_lengths, results.injected_burst_count);
grid on;
xlabel("Injected burst length");
ylabel("Count");
title(sprintf("PAM-%d, F(D)=(1+D)^%d", M, pr_order));
nexttile;
bar(1:numel(results.output_burst_count), results.output_burst_count);
grid on;
xlabel("Output error burst length");
ylabel("Count");
title("Output symbol-error burst count");
end
function printBurstCountReport(results, M, pr_order)
%PRINTBURSTCOUNTREPORT Paste-friendly burst-count report.
max_len = max([numel(results.input_burst_count), numel(results.output_burst_count)]);
lengths = 1:max_len;
input_counts = padCountVector(results.input_burst_count, max_len);
output_counts = padCountVector(results.output_burst_count, max_len);
fprintf('BURST_COUNT_REPORT pr_order=%d M=%d\n', pr_order, M);
fprintf('length : %s\n', formatIntegerVector(lengths));
fprintf('before_cnt : %s\n', formatIntegerVector(input_counts));
fprintf('after_cnt : %s\n', formatIntegerVector(output_counts));
fprintf('delta_cnt : %s\n', formatIntegerVector(output_counts - input_counts));
fprintf('before_total=%d before_max=%d after_total=%d after_max=%d\n', ...
results.input_burst_total, results.input_max_burst_length, ...
results.output_burst_total, results.output_max_burst_length);
end
function x = padCountVector(x, target_length)
%PADCOUNTVECTOR Pad count vector with trailing zeros.
x = reshape(x, 1, []);
if numel(x) < target_length
x(end+1:target_length) = 0;
end
end
function s = formatIntegerVector(x)
%FORMATINTEGERVECTOR Format vector as comma-separated integers.
x = reshape(x, 1, []);
if isempty(x)
s = '';
return
end
parts = strings(size(x));
for k = 1:numel(x)
parts(k) = string(sprintf('%d', x(k)));
end
s = char(strjoin(parts, ', '));
end
function plotBurstSummary(burst_results, pr_orders, M_list)
%PLOTBURSTSUMMARY Show BER and output-burst metrics for every M/order pair.
ber_map = nan(numel(pr_orders), numel(M_list));
symbol_error_map = nan(size(ber_map));
max_burst_map = nan(size(ber_map));
for row = 1:numel(pr_orders)
pr_order = pr_orders(row);
for col = 1:numel(M_list)
M = M_list(col);
results = burst_results(pr_order, M).results;
ber_map(row, col) = results.ber;
symbol_error_map(row, col) = results.symbol_errors;
max_burst_map(row, col) = results.output_max_burst_length;
end
end
figure(500);
clf;
tiledlayout(1, 3);
nexttile;
imagesc(M_list, pr_orders, ber_map);
colorbar;
xlabel("PAM order M");
ylabel("PR order n");
title("BER");
nexttile;
imagesc(M_list, pr_orders, symbol_error_map);
colorbar;
xlabel("PAM order M");
ylabel("PR order n");
title("Symbol errors");
nexttile;
imagesc(M_list, pr_orders, max_burst_map);
colorbar;
xlabel("PAM order M");
ylabel("PR order n");
title("Max output burst");
end
function [u, y_amp, pr_alphabet, y_idx, d_ref] = emulatePrecoder(d_tx, M, init_state, pr_order)
%EMULATEPRECODER PR precoder emulation by reference relabeling.
% The available sequence d_tx is interpreted as d_pre.
% d_ref(k) = sum_i h(i+1)*d_tx(k-i) mod M
% The waveform itself is generated by the ordinary precode/encode blocks.
if nargin < 4 || isempty(pr_order)
pr_order = 1;
end
if nargin < 3
init_state = [];
end
d_ref = precoderEmulationReference(d_tx, M, init_state, pr_order);
u = precode(d_ref, M, init_state, pr_order);
[y_amp, pr_alphabet, y_idx] = encode(u, M, init_state, pr_order);
assert(isequal(u, d_tx), ...
'Emulation relabeling did not reproduce the requested precoder output.');
end
function d_ref = precoderEmulationReference(d_tx, M, init_state, pr_order)
%PRECODEREMULATIONREFERENCE Reference implied by interpreting d_tx as u.
% For h(D)=(1+D)^n, d_ref(k)=mod(sum_i h(i+1) d_tx(k-i), M).
if nargin < 4 || isempty(pr_order)
pr_order = 1;
end
if nargin < 3
init_state = [];
end
validateSymbolVector(d_tx, M, 'd_tx');
pr_order = validatePrOrder(pr_order);
state = validateInitState(init_state, M, pr_order);
original_size = size(d_tx);
d_work = double(d_tx(:).');
h = getPartialResponseTaps(pr_order);
d_ref = zeros(size(d_work));
for k = 1:numel(d_work)
pr_state = [d_work(k) state];
d_ref(k) = mod(sum(h .* pr_state), M);
state = pr_state(1:end-1);
end
d_ref = reshape(d_ref, original_size);
end
function s = symbolIndexToAmplitude(a, M)
%SYMBOLINDEXTOAMPLITUDE Map natural PAM indices to amplitudes.
validateSymbolVector(a, M, 'a');
pam_levels = getPamLevels(M);
s = reshape(pam_levels(double(a(:)) + 1), size(a));
end
function a = amplitudeToSymbolIndex(s, M)
%AMPLITUDETOSYMBOLINDEX Map exact PAM amplitudes to natural indices.
validatePamOrder(M);
validateRealVector(s, 's');
original_size = size(s);
s_work = double(s(:).');
pam_levels = getPamLevels(M);
a = zeros(size(s_work));
for k = 1:numel(s_work)
match = find(pam_levels == s_work(k), 1);
if isempty(match)
error('amplitudeToSymbolIndex:InvalidAmplitude', ...
'All amplitudes must be valid PAM-%d levels.', M);
end
a(k) = match - 1;
end
a = reshape(a, original_size);
end
function runDuobinaryPartialResponseTests()
%RUNDUOBINARYPARTIALRESPONSETESTS Minimal self-checks.
fprintf('\nRunning duobinary PR self-tests...\n');
testNoiselessReconstruction(2, 1);
testNoiselessReconstruction(4, 1);
testNoiselessReconstruction(8, 1);
testNoiselessReconstruction(4, 2);
testNoiselessReconstruction(8, 3);
testFrameworkNoiselessBer();
testFirstOrderPrecoderEmulation();
testAlphabetSize();
testNonzeroInitState();
fprintf('All duobinary PR self-tests passed.\n');
end
function testFrameworkNoiselessBer()
for M = [2 4 8]
for pr_order = [1 2 3]
[txBits, txSymbols, mapper, a] = generatePamSymbols(M, 1024, 9);
u = precode(a, M, [], pr_order);
y_amp = encode(u, M, [], pr_order);
a_hat = memorylessdecode(y_amp, M, [], pr_order);
rxBits = demapSymbolIndices(a_hat, txSymbols, mapper, M);
[~, errors, ber] = calc_ber(rxBits.signal, txBits.signal);
assert(errors == 0 && ber == 0, ...
'Framework BER check failed for M=%d, order=%d.', M, pr_order);
end
end
fprintf('PASS: framework Signalgenerator/PAMmapper/calc_ber noiseless BER\n');
end
function testFirstOrderPrecoderEmulation()
for M = [2 4 8]
for pr_order = [1 2 3]
init_state = mod((1:pr_order) + M - 2, M);
[~, txSymbols, mapper, d_tx] = generatePamSymbols(M, 1024, 9);
[u, y_amp, ~, ~, d_ref] = emulatePrecoder(d_tx, M, init_state, pr_order);
assert(isequal(u, d_tx), ...
'Emulated precoder output mismatch for M=%d, order=%d.', M, pr_order);
d_hat = memorylessdecode(y_amp, M, init_state, pr_order);
assert(isequal(d_hat, d_ref), ...
'Emulated noiseless modulo decoding failed for M=%d, order=%d.', M, pr_order);
refBits = demapSymbolIndices(d_ref, txSymbols, mapper, M);
rxBits = demapSymbolIndices(d_hat, txSymbols, mapper, M);
[~, errors, ber] = calc_ber(rxBits.signal, refBits.signal);
assert(errors == 0 && ber == 0, ...
'Emulated framework BER failed for M=%d, order=%d.', M, pr_order);
end
end
fprintf('PASS: precoder emulation and relabeled reference\n');
end
function testNoiselessReconstruction(M, pr_order)
a = randi([0 M-1], 1000, 1);
u = precode(a, M, [], pr_order);
y_amp = encode(u, M, [], pr_order);
a_hat = memorylessdecode(y_amp, M, [], pr_order);
assert(isequal(a_hat, a), ...
'Noiseless reconstruction failed for M=%d, order=%d.', M, pr_order);
fprintf('PASS: noiseless reconstruction for M=%d, order=%d\n', M, pr_order);
end
function testAlphabetSize()
for M = [2 4 8]
for pr_order = [1 2 3]
u = precode(0:M-1, M, [], pr_order);
[~, pr_alphabet] = encode(u, M, [], pr_order);
expected_size = (M-1)*2^pr_order + 1;
assert(numel(pr_alphabet) == expected_size, ...
'PR alphabet size failed for M=%d, order=%d.', M, pr_order);
end
end
fprintf('PASS: PR alphabet has (M-1)*2^order+1 distinct levels\n');
end
function testNonzeroInitState()
for M = [2 4 8]
for pr_order = [1 2 3]
init_state = mod((1:pr_order) + M - 2, M);
a = randi([0 M-1], 1, 1000);
u = precode(a, M, init_state, pr_order);
y_amp = encode(u, M, init_state, pr_order);
a_hat = memorylessdecode(y_amp, M, init_state, pr_order);
assert(isequal(a_hat, a), ...
'Nonzero init_state failed for M=%d, order=%d.', M, pr_order);
end
end
fprintf('PASS: random-state test with nonzero init_state\n');
end
function [pr_alphabet, idx_for_alphabet] = getDuobinaryAlphabet(M, pr_order)
%GETDUOBINARYALPHABET Build amplitude/index lookup tables from PAM sums.
validatePamOrder(M);
pr_order = validatePrOrder(pr_order);
h = getPartialResponseTaps(pr_order);
pam_levels = getPamLevels(M);
n_state = pr_order + 1;
n_combinations = M^n_state;
all_amp = zeros(1, n_combinations);
idx_to_amp = nan(1, (M-1)*sum(h) + 1);
for c = 0:n_combinations-1
state = zeros(1, n_state);
tmp = c;
for k = 1:n_state
state(k) = mod(tmp, M);
tmp = floor(tmp / M);
end
y_idx = sum(h .* state);
y_amp = sum(h .* pam_levels(state + 1));
all_amp(c + 1) = y_amp;
idx_to_amp(y_idx + 1) = y_amp;
end
pr_alphabet = unique(all_amp);
idx_for_alphabet = zeros(size(pr_alphabet));
for k = 1:numel(pr_alphabet)
idx_for_alphabet(k) = find(idx_to_amp == pr_alphabet(k), 1) - 1;
end
end
function h = getPartialResponseTaps(pr_order)
%GETPARTIALRESPONSETAPS Coefficients of F(D)=(1+D)^n.
pr_order = validatePrOrder(pr_order);
h = zeros(1, pr_order + 1);
for k = 0:pr_order
h(k + 1) = nchoosek(pr_order, k);
end
end
function validatePamOrder(M)
if ~(isscalar(M) && isnumeric(M) && any(M == [2 4 8]))
error('validatePamOrder:InvalidM', 'M must be one of [2 4 8].');
end
end
function validateSymbolVector(a, M, name)
validatePamOrder(M);
if ~(isnumeric(a) && isreal(a) && isvector(a))
error('validateSymbolVector:InvalidVector', '%s must be a real numeric vector.', name);
end
if any(a(:) ~= floor(a(:))) || any(a(:) < 0) || any(a(:) > M-1)
error('validateSymbolVector:InvalidSymbols', ...
'%s must contain integers in [0, M-1].', name);
end
end
function validateRealVector(x, name)
if ~(isnumeric(x) && isreal(x) && isvector(x))
error('validateRealVector:InvalidVector', '%s must be a real numeric vector.', name);
end
end
function pr_order = validatePrOrder(pr_order)
if ~(isscalar(pr_order) && isnumeric(pr_order) && isreal(pr_order) && ...
pr_order == floor(pr_order) && pr_order >= 1)
error('validatePrOrder:InvalidOrder', ...
'pr_order must be a positive integer.');
end
pr_order = double(pr_order);
end
function init_state = validateInitState(init_state, M, pr_order)
validatePamOrder(M);
pr_order = validatePrOrder(pr_order);
if isempty(init_state)
init_state = zeros(1, pr_order);
end
if ~(isnumeric(init_state) && isreal(init_state) && isvector(init_state) && ...
numel(init_state) == pr_order)
error('validateInitState:InvalidState', ...
'init_state must contain pr_order integers in [0, M-1].');
end
if any(init_state(:) ~= floor(init_state(:))) || ...
any(init_state(:) < 0) || any(init_state(:) > M-1)
error('validateInitState:InvalidState', ...
'init_state must contain pr_order integers in [0, M-1].');
end
init_state = double(reshape(init_state, 1, []));
end

0
Theory/Mapping_Coding/partial_response/partial_response_filter.m → Theory/Mapping_Coding/partial_response/plot_partial_response_filter.m
View File

4
Theory/Mapping_Coding/partial_response/test_db_minimal_example.m
View File

@@ -1,5 +1,5 @@
M = 4;
M = 8;
apply_precode_at_tx = 1;
bits = Signalgenerator( ...
@@ -19,9 +19,7 @@ disp(['Tx Sequenz: -- RMS:',sprintf('%.1f',rms(symbols_tx.signal)),' - - Levels
unique(symbols_tx.signal)
disp('- - - - - - - - - -');
symbols_tx.signal = awgn(symbols_tx.signal,20,"measured",1);
show2Dconstellation(symbols_tx,symbols_tx,"displayname",'VNLE Out','fignum',2241);
if apply_precode_at_tx
% Entschiedene Symbole codieren: d_DB(n) = d(n) + d(n-1) (im Fall von PAM4 7 level [0 1 2 3 4 5 6])

19
Theory/Mapping_Coding/signalgenerator_pam_ber_minimal_example.m
View File

@@ -12,29 +12,16 @@ fprintf("PAM | bits checked | errors | BER\n");
fprintf("----+--------------+--------+-----\n");
for M = orders
if M == 6
bitsPerBlock = 5;
else
bitsPerBlock = log2(M);
end
bitSignal = Signalgenerator( ...
"form", signalform.prms, ...
"dimension", bitsPerBlock, ...
"M", M, ...
"order", userOrder).process();
if M == 6
txBits = reshape(bitSignal.signal.', [], 1);
txBits = txBits(1:end - mod(numel(txBits), bitsPerBlock));
else
txBits = bitSignal.signal;
end
mapper = PAMmapper(M, 0);
txSymbols = mapper.map(txBits);
txSymbols = mapper.map(bitSignal);
rxBits = mapper.demap(txSymbols);
[checkedBits, errors, ber] = calc_ber(rxBits, txBits);
[checkedBits, errors, ber] = calc_ber(rxBits.signal, bitSignal.signal);
fprintf("%3d | %12d | %6d | %.1e\n", M, checkedBits, errors, ber);
end