From 456fd02186ba69532054a2b3118e4fa314e2f654 Mon Sep 17 00:00:00 2001 From: Silas Oettinghaus Date: Thu, 23 Apr 2026 16:16:58 +0200 Subject: [PATCH] Ordnerstruktur --- Classes/00_signals/Signal.m | 154 ++-- Classes/01_transmit/Signalgenerator.m | 64 +- Functions/EQ_visuals/showErrorBurstCount.m | 36 +- Functions/Metrics/count_error_bursts.m | 6 +- Tests/00_signals/Signal_test.m | 47 + Tests/01_transmit/Signalgenerator_test.m | 29 + .../partial_response/duobinary_back2back.m | 18 + .../duobinary_minimal_example.m | 26 +- .../pam_6_differential_code_understand.m | 73 -- .../partial_response_precoder.m | 835 ++++++++++++++++++ ...ilter.m => plot_partial_response_filter.m} | 0 .../test_db_minimal_example.m | 4 +- .../signalgenerator_pam_ber_minimal_example.m | 19 +- 13 files changed, 1117 insertions(+), 194 deletions(-) create mode 100644 Theory/Mapping_Coding/partial_response/duobinary_back2back.m delete mode 100644 Theory/Mapping_Coding/partial_response/pam_6_differential_code_understand.m create mode 100644 Theory/Mapping_Coding/partial_response/partial_response_precoder.m rename Theory/Mapping_Coding/partial_response/{partial_response_filter.m => plot_partial_response_filter.m} (100%) diff --git a/Classes/00_signals/Signal.m b/Classes/00_signals/Signal.m index 407691c..3579fa5 100644 --- a/Classes/00_signals/Signal.m +++ b/Classes/00_signals/Signal.m @@ -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); diff --git a/Classes/01_transmit/Signalgenerator.m b/Classes/01_transmit/Signalgenerator.m index 3062a07..4c9a5ad 100644 --- a/Classes/01_transmit/Signalgenerator.m +++ b/Classes/01_transmit/Signalgenerator.m @@ -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 diff --git a/Functions/EQ_visuals/showErrorBurstCount.m b/Functions/EQ_visuals/showErrorBurstCount.m index cb9e016..649b88b 100644 --- a/Functions/EQ_visuals/showErrorBurstCount.m +++ b/Functions/EQ_visuals/showErrorBurstCount.m @@ -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 \ No newline at end of file +end diff --git a/Functions/Metrics/count_error_bursts.m b/Functions/Metrics/count_error_bursts.m index 47d5439..6fd38b4 100644 --- a/Functions/Metrics/count_error_bursts.m +++ b/Functions/Metrics/count_error_bursts.m @@ -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); diff --git a/Tests/00_signals/Signal_test.m b/Tests/00_signals/Signal_test.m index 829ee1c..6bb3dda 100644 --- a/Tests/00_signals/Signal_test.m +++ b/Tests/00_signals/Signal_test.m @@ -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]); diff --git a/Tests/01_transmit/Signalgenerator_test.m b/Tests/01_transmit/Signalgenerator_test.m index d4a6337..715aa32 100644 --- a/Tests/01_transmit/Signalgenerator_test.m +++ b/Tests/01_transmit/Signalgenerator_test.m @@ -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, ... diff --git a/Theory/Mapping_Coding/partial_response/duobinary_back2back.m b/Theory/Mapping_Coding/partial_response/duobinary_back2back.m new file mode 100644 index 0000000..d2b6e33 --- /dev/null +++ b/Theory/Mapping_Coding/partial_response/duobinary_back2back.m @@ -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 \ No newline at end of file diff --git a/Theory/Mapping_Coding/partial_response/duobinary_minimal_example.m b/Theory/Mapping_Coding/partial_response/duobinary_minimal_example.m index b8f0703..71e26a3 100644 --- a/Theory/Mapping_Coding/partial_response/duobinary_minimal_example.m +++ b/Theory/Mapping_Coding/partial_response/duobinary_minimal_example.m @@ -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 \ No newline at end of file +legend diff --git a/Theory/Mapping_Coding/partial_response/pam_6_differential_code_understand.m b/Theory/Mapping_Coding/partial_response/pam_6_differential_code_understand.m deleted file mode 100644 index c019c12..0000000 --- a/Theory/Mapping_Coding/partial_response/pam_6_differential_code_understand.m +++ /dev/null @@ -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)'; \ No newline at end of file diff --git a/Theory/Mapping_Coding/partial_response/partial_response_precoder.m b/Theory/Mapping_Coding/partial_response/partial_response_precoder.m new file mode 100644 index 0000000..340cc8a --- /dev/null +++ b/Theory/Mapping_Coding/partial_response/partial_response_precoder.m @@ -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 + burst_label_plan(end + 1, 1) = burst_length; %#ok + 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 + err_label = [err_label; repmat(burst_label_plan(b), numel(current_pos), 1)]; %#ok + 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 diff --git a/Theory/Mapping_Coding/partial_response/partial_response_filter.m b/Theory/Mapping_Coding/partial_response/plot_partial_response_filter.m similarity index 100% rename from Theory/Mapping_Coding/partial_response/partial_response_filter.m rename to Theory/Mapping_Coding/partial_response/plot_partial_response_filter.m diff --git a/Theory/Mapping_Coding/partial_response/test_db_minimal_example.m b/Theory/Mapping_Coding/partial_response/test_db_minimal_example.m index e0ef0af..cfeb77d 100644 --- a/Theory/Mapping_Coding/partial_response/test_db_minimal_example.m +++ b/Theory/Mapping_Coding/partial_response/test_db_minimal_example.m @@ -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]) diff --git a/Theory/Mapping_Coding/signalgenerator_pam_ber_minimal_example.m b/Theory/Mapping_Coding/signalgenerator_pam_ber_minimal_example.m index 8868ad9..559a844 100644 --- a/Theory/Mapping_Coding/signalgenerator_pam_ber_minimal_example.m +++ b/Theory/Mapping_Coding/signalgenerator_pam_ber_minimal_example.m @@ -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