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

New MPI mitigation schemes // Duobinary // Start of FTN schemes

Browse Source
This commit is contained in:
Silas Oettinghaus
2024-09-06 08:10:04 +02:00
parent bb228ae2bd
commit 03bfd70470
30 changed files with 1598 additions and 560 deletions
Download Patch File Download Diff File Expand all files Collapse all files

114
Classes/00_signals/Signal.m
View File

@@ -237,12 +237,16 @@ classdef Signal
arguments
obj Signal
options.fs_in double
options.fs_in double = obj.fs
options.fs_out double
options.n double = 10;
options.beta double = 5;
end
if options.fs_in ~= obj.fs
warning('The signals fs is different from the given fs_in while it should be the same.');
end
obj.signal = resample(obj.signal,options.fs_out,options.fs_in,options.n,options.beta);
desc = ['resample signal from ', num2str(options.fs_in*1e-9), ' GHz to ', num2str(options.fs_out*1e-9), ' GHz' ];
@@ -265,7 +269,7 @@ classdef Signal
% spectrum_plot(obj.signal,options.fsamp,options.figurename,options.displayname);
N = 2^(nextpow2(length(obj.signal))-8);
[p_lin,w] = pwelch(obj.signal,hanning(N),N/2,N,obj.fs,"centered","power","mean");
p_dbm = 10*log10(p_lin)+30; %dB to dBm in case of "power"
@@ -279,7 +283,7 @@ classdef Signal
edgetick = 2^(nextpow2(obj.fs*1e-9));
% xticks([-edgetick:16:edgetick]);
xlim([-244, 244])
ylim([-120,-0]);
ylim([-120,10]);
yticks([-200:10:10]);
legend
@@ -301,7 +305,7 @@ classdef Signal
pow = pow / 50;
end
pow = 10*log10(pow)+30; %dbm
case power_notation.mW
pow = pow .* 1e3; %mW
@@ -378,14 +382,23 @@ classdef Signal
end
%%
function [obj,delay_n] = delay(obj,options)
function [obj] = delay(obj,delay,options)
arguments
obj Signal
options.delay_samples = 0
delay double = 0
options.mode delay_mode = delay_mode.samples
end
obj.signal=delayseq(obj.signal,options.delay_samples);
if options.mode == delay_mode.samples
obj.signal=delayseq(obj.signal,delay);
elseif options.mode == delay_mode.time
obj.signal=delayseq(obj.signal,delay,obj.fs);
end
end
@@ -415,7 +428,7 @@ classdef Signal
end
% delay by lagging samples
obj = obj.delay("delay_samples",-D);
obj = obj.delay(-D,'mode','samples');
cuts = obj.length-(options.reference.length*q);
@@ -491,7 +504,7 @@ classdef Signal
ppeak(i) = maxA - (difference/histpoints*loc(i));
end
if isa(obj,'Opticalsignal')
er=10*log10(ppeak(1)/ppeak(end));
@@ -633,31 +646,7 @@ classdef Signal
hold on
xline(posxall)
hist_interest = plot_data(:,posxall);
hist_interest_smoth = smooth(hist_interest,20);
a = scatter(hist_interest_smoth+posxall,1:length(hist_interest_smoth),4,'.','MarkerEdgeColor','red');
[pk,loc] = findpeaks(hist_interest_smoth,"MinPeakDistance",40,"NPeaks",M,"MinPeakHeight",30);
scatter(posxall,loc,'red','Marker','x','LineWidth',2);
yline(loc,'Color','red','LineWidth',1,'LineStyle',':');
for i = 1:numel(loc)
ppeak(i) = maxA - (difference/histpoints*loc(i));
end
oma = false;
if isa(obj,'Opticalsignal')
er=10*log10(ppeak(1)/ppeak(end));
elseif isa(obj,'Electricalsignal')
if mean([ppeak(1),ppeak(end)]) < 1e-2
oma = true;
er=max(ppeak)-min(ppeak);
else
er=10*log10(ppeak(1)/ppeak(end));
end
else
er=10*log10(ppeak(1)/ppeak(end));
end
% Define properties
boxPosition = [0.15 0.86 0.2 0.05]; % Position for the first box [x y width height]
@@ -687,24 +676,53 @@ classdef Signal
'FontWeight', boxFontWeight, ...
'HorizontalAlignment', 'center');
% Adjust position for the third box (slightly to the right)
boxPosition = [0.59 0.86 0.2 0.05]; % Adjusted position
try
hist_interest = plot_data(:,posxall);
hist_interest_smoth = smooth(hist_interest,20);
a = scatter(hist_interest_smoth+posxall,1:length(hist_interest_smoth),4,'.','MarkerEdgeColor','red');
% Create third annotation box for Vmax
if ~oma
thirdboxstring = ['ER (db):',num2str(er),' dB'];
else
thirdboxstring = ['OMA outer:',num2str(er),' V'];
end
[pk,loc] = findpeaks(hist_interest_smoth,"MinPeakDistance",40,"NPeaks",M,"MinPeakHeight",30);
scatter(posxall,loc,'red','Marker','x','LineWidth',2);
yline(loc,'Color','red','LineWidth',1,'LineStyle',':');
annotation('textbox', boxPosition, ...
'String',thirdboxstring , ...
'BackgroundColor', boxColor, ...
'EdgeColor', boxEdgeColor, ...
'LineStyle', boxLineStyle, ...
'FontWeight', boxFontWeight, ...
'HorizontalAlignment', 'center');
for i = 1:numel(loc)
ppeak(i) = maxA - (difference/histpoints*loc(i));
end
oma = false;
if isa(obj,'Opticalsignal')
er=10*log10(ppeak(1)/ppeak(end));
elseif isa(obj,'Electricalsignal')
if mean([ppeak(1),ppeak(end)]) < 1e-2
oma = true;
er=max(ppeak)-min(ppeak);
else
er=10*log10(ppeak(1)/ppeak(end));
end
else
er=10*log10(ppeak(1)/ppeak(end));
end
% Adjust position for the third box (slightly to the right)
boxPosition = [0.59 0.86 0.2 0.05]; % Adjusted position
% Create third annotation box for Vmax
if ~oma
thirdboxstring = ['ER (db):',num2str(er),' dB'];
else
thirdboxstring = ['OMA outer:',num2str(er),' V'];
end
annotation('textbox', boxPosition, ...
'String',thirdboxstring , ...
'BackgroundColor', boxColor, ...
'EdgeColor', boxEdgeColor, ...
'LineStyle', boxLineStyle, ...
'FontWeight', boxFontWeight, ...
'HorizontalAlignment', 'center');
end
yticks(linspace(0,histpoints,16));
y_tickstring = sprintfc('%.2f', y_tickstring);

406
Classes/01_transmit/ChannelFreqResp.m Normal file
View File

@@ -0,0 +1,406 @@
classdef ChannelFreqResp < handle
% Linear pre-compensation of Channel effects. This module has two modes:
% A) In acquire mode it sends a real OFDM singal with all subcarriers assigned in
% order to acquire later the system's frequency response after the
% signal passed the system...
% ChannelFreqResp.SendOFDM
% B) In the "not acquire" mode it loads the frequency response generated by the freq_res module,
% calculates the inverse frequency response with the right length matched
% with the signal length and distorts the signal before tranmission
% give channel indices for measure mode
%TYPICAL FLOW:
% referencesignal = ChannelFreqResp.buildOFDM()
% referencesignal -> SYSTEM -> measuredsignal
% ChannelFreqResp.estimate(measuredsignal)
properties(Access=public)
Nacq
Navg
Ncp
f_ref %fs of transmitted dmt sig
f_observe %fs of the observed signal
symlen
seqlen
refsig
H
H_all
H_apply
Nfft
df
faxis
end
methods (Access=public)
function obj = ChannelFreqResp(options)
%NAME Construct an instance of this class
% Detailed explanation goes here
arguments
options.Nacq = 4096;
options.Navg = 30;
options.Ncp = 100;
options.f_ref;
options.f_observe;
end
%
fn = fieldnames(options);
for n = 1:numel(fn)
try
obj.(fn{n}) = options.(fn{n});
end
end
obj.Nfft = 2*obj.Nacq + 1 ;
obj.df = obj.f_ref/obj.Nfft ; % calculate frequency grid
obj.faxis = (0:obj.Nfft-1)*obj.df ;
end
function output = buildOFDM(obj)
% generate OFDM or DMT signal that is send over the channel
obj.symlen = (obj.Nacq*2+obj.Ncp+1);
obj.seqlen = obj.symlen*obj.Navg;
obj.refsig = obj.genRNDDMT(15)';
obj.refsig = Informationsignal(obj.refsig,"fs",obj.f_ref);
output = obj.refsig;
end
function output = estimate(obj, data_in, options)
arguments
obj
data_in
options.save logical = false
options.savePath = "";
options.fileName = "";
end
% use the transmitted DMT and the received DMT to get the transfer function
obj.Nfft = 2*obj.Nacq + 1 ;
%resample from fadc to fdac (fdac => fs of reference signal)
data_in = data_in.resample("fs_out",obj.f_ref);
%est from dmt
obj.H_all = obj.estHfromDMT(data_in.signal,obj.refsig.signal);
obj.H = mean(obj.H_all,1);
output = obj.H;
if options.save
obj.save('fileName',options.fileName,'savePath',options.savePath);
end
end
function Target = precomp(obj,Target,options)
% apply the freq response
arguments
obj
Target
options.fileName = ''
options.loadPath = ''
options.maxampdb
end
if isempty(obj.H)
obj.load("fileName",options.fileName,"loadPath",options.loadPath);
end
H_inv = 1./obj.H;
nH = find(~isnan(H_inv),1,'last'); %last value that is not nan
H_inv(nH:end)=H_inv(nH); %replace everything after first NaN with the last non-NaN value
fstarget = Target.fs;
% Build new frequencie axis (with current fs)
fnew = linspace(0,fstarget/2,length(Target)/2+1);
fnew = fnew(2:end-1);
% Old frequency axis (should be much coarser)
idx_old = find((obj.faxis > 0) .* (obj.faxis < fstarget/2)); %positions of all Frequencies smaller than fs/2
int_fold = obj.faxis(idx_old); %old frequencies from 0 to fs/2
idx_new = find(fnew <= int_fold(end)); %positions of all Frequencies smaller than fs/2
int_fnew = fnew(idx_new);%new frequencies from 0 to fs/2
% interpolate the frequency response that had a coarse frequency resolution (e.g. 256 bins) to the current frequency resolution (e.g. 21843 bins)
iH = interp1(int_fold, real(H_inv(idx_old)) ,fnew, 'linear') + 1i*interp1(int_fold, imag(H_inv(idx_old)) ,fnew, 'linear');
% set all NaN values to the fist/ last non-NaN value
nH = find(~isnan(iH),1,'first');
iH(1:nH)=iH(nH);
nH = find(~isnan(iH),1,'last');
iH(nH:end) =iH(nH);
%smoothing takes time and sometimes the result looks odd,
%however the performance is most of the time better
smoothing = 1;
if smoothing
iH = smooth(fnew,iH,0.1,'loess')';
end
% multiply the complex frequency responce with the phase
% angle (-pi,pi) at lowest frequency
iH = iH.*exp(-1j*angle(iH(1))); % to be checked (<- not from silas, so what needs to be checked?)
% Phase difference between lowest and hiughest frequency
% component -> but what is this for? dPhase is not used...
noncausal = 0;
if noncausal
dPhase = angle(iH(end))-angle(iH(1));
end
% normalize complex freq. resp. by magnitude at the first
% five frequencies -> should be the vaue at f=0=DC component?
iH = iH./mean(abs(iH(1:100))); %why 1:5??
% set maximum amplification
% set als values higher than hmax to hmax and keep the
% phase information by multiplication with respective
% corresponding phase angles
maxamp_lin = 10^(options.maxampdb/20);
iH(abs(iH)>maxamp_lin) = maxamp_lin.*exp(1j*angle(iH(abs(iH)>maxamp_lin)));
% it could be helpful to clip at linear 1 (=to keep comp from attenuating the signal)
% iH(abs(iH)<1) = 1.*exp(1j*angle(iH(abs(iH)<1)));
if 0
figure(7);hold on;plot(fnew,20*log10(abs(iH)))
end
H_inv = [iH(1) iH fliplr(conj(iH)) conj(iH(1))];
obj.H_apply = H_inv;
Target.signal = real((ifft( ( fft(real( Target.signal )) .* H_inv' ) )));
end
function plot(obj)
figure(55551);
clf;
Havg = obj.H;
%1)
subplot(4,1,1);hold all;box on;title('Magnitude Freq. Response');
plot(obj.faxis/1e9, 20*log10(abs(obj.H_all)),'linewidth',0.1,'LineStyle','-','Color','#808080') ;
xlim([0.2 .5*max(obj.faxis)*1e-9]);
plot(obj.faxis/1e9, 20*log10(abs(Havg)),'LineWidth',2);
grid on;
%2)
subplot(4,1,2); hold all; box on; title('Phase Freq. Response');
plot(obj.faxis/1e9, (angle(obj.H_all)),'linewidth',0.1,'LineStyle','-','Color','#808080') ;
plot(obj.faxis/1e9, (angle(Havg)),'LineWidth',2) ;
xlim([0.2 .5*max(obj.faxis)*1e-9]);
grid on;
%normalize / remove attenuation
Havg = Havg./mean(Havg(2:10));
%3)
subplot(4,1,3); hold all; box on; title('Inverse Magnitude Freq. Response');
plot(obj.faxis/1e9, 20*log10(abs(1./Havg)),"LineWidth",2,"Color",[0.3467 0.5360 0.6907]) ;
xlim([0.2 .5*max(obj.faxis)*1e-9]); grid on;
ylim([-1 15]);
hold on;
yline(3,'LineWidth',2,'LineStyle','--');
%4)
subplot(4,1,4); hold all; box on; title('Inverse Phase Freq. Response');
plot(obj.faxis/1e9, (angle(1./Havg)),"LineWidth",2,"Color",[0.3467 0.5360 0.6907]) ;
xlim([0.2 .5*max(obj.faxis)*1e-9]); grid on;
end
function save(obj,options)
arguments
obj
options.fileName = ''
options.savePath = ''
end
% Check if the fileName was provided
if isempty(char(options.fileName))
% If no file name is provided, prompt the user to enter a file name
[options.fileName, filePath] = uiputfile('*.mat', 'Save As');
% If the user cancels the dialog, fileName and filePath will be 0
if isequal(options.fileName, 0) || isequal(filePath, 0)
disp('Save operation cancelled.');
return;
end
% Update the savePath with the directory chosen by the user
options.savePath = filePath;
end
% Check if the savePath was provided
if isempty(char(options.savePath))
% If no path is provided, open a UI window to select the path
options.savePath = uigetdir('', 'Select a folder to save the file');
% If the user cancels the dialog, savePath will be 0
if options.savePath == 0
disp('Save operation cancelled.');
return;
end
else
% If a path is provided, validate it
if ~isfolder(options.savePath)
error('The specified path does not exist.');
end
end
% Construct the full file path
fullFileName = fullfile(options.savePath, options.fileName);
% Save the data to the specified file
save(fullFileName, 'obj');
fprintf('Frequency response information successfully saved to %s\n', fullFileName);
end
function data = load(obj, options)
% Function to load data from a specified file and path.
arguments
obj
options.fileName = ''
options.loadPath = ''
end
% Check if the fileName was provided
if isempty(char(options.fileName))
% If no file name is provided, open a UI window to select the file
[options.fileName, options.loadPath] = uigetfile('*.mat', 'Select a file to load');
% If the user cancels the dialog, fileName and loadPath will be 0
if isequal(options.fileName, 0) || isequal(options.loadPath, 0)
disp('Load operation cancelled.');
data = [];
return;
end
end
% If loadPath is not provided or empty, use the current folder
if isempty(char(options.loadPath))
options.loadPath = pwd;
else
% Validate the path if provided
if ~isfolder(options.loadPath)
error('The specified path does not exist.');
end
end
% Construct the full file path
fullFileName = fullfile(options.loadPath, options.fileName);
% Check if the file exists
if ~isfile(fullFileName)
fullFileName = fullfile(options.loadPath, [char(options.fileName),'.mat']);
if ~isfile(fullFileName)
error('The specified file does not exist.');
end
end
% Load the data from the specified file
loadedData = load(fullFileName);
% Replace whole obj here.. is this save or unsave?!
fn = fieldnames(loadedData.obj);
for n = 1:numel(fn)
try
obj.(fn{n}) = loadedData.obj.(fn{n});
end
end
fprintf('Frequency response information successfully loaded from %s\n', fullFileName);
end
end
methods (Access=private)
% Cant be seen from outside! So put all your functions here that can/
% shall not be called from outside
function rOFDM = genRNDDMT(obj,randkey)
rOFDM = NaN(1,(2*obj.Nacq+obj.Ncp+1)*obj.Navg);
s = RandStream('mt19937ar','Seed',randkey,'NormalTransform','Polar');
ref = 2 * round(rand(s,obj.Nacq, obj.Navg)) - 1 ; % Navg Random BPSK sequences
ofdm = [ones(1, obj.Navg) ; ref ; conj(ref(end:-1:1,:)) ] ; % DMT
ofdm = ifft(ofdm) ; % Navg real OFDM sequences
ofdm = [ofdm(end-obj.Ncp+1 : end, :) ; ofdm] ; % Add cyclic prefix
rOFDM = reshape(ofdm, 1, size(ofdm,1)*size(ofdm,2)) ;
rOFDM = rOFDM./max(abs(rOFDM(:)));
end
function [rH] = estHfromDMT(obj, data_in, ref_in)
%! Dont (circ)shift the signal here as this would remove the phase information!
%tested with a butterworth filter this exactly reconstructs the
%phase and the magnitude. However, the option is here
estimatephase = 1;
if ~estimatephase
% 0. cross-correlate the received signal with its reference to extract the periods
corr = abs(ifft( fft(data_in(1:length(ref_in))).* conj(fft(ref_in)) )) ;
% find max
[~, peak] = max(corr) ;
peak=max(1,peak-1);
data_in = circshift(data_in,-peak) ;
%Y = data_in(peak:peak+(Nfft+obj.Ncp)*obj.Navg-1) ;
end
% 1. Reshape signal to a matrix to support noise averaging
Nfft = 2*obj.Nacq + 1 ;
Y = reshape(data_in, Nfft+obj.Ncp, obj.Navg).' ;
X = reshape(ref_in, Nfft + obj.Ncp, obj.Navg).' ;
% 2. Remove cyclic prefix and apply FFT transformation
Y = fft(Y(:, obj.Ncp+1 : obj.Ncp + Nfft), [], 2) ;
X = fft(X(:, obj.Ncp+1 : obj.Ncp + Nfft), [], 2) ;
% 4. Caclulate the frequency response using H = Y/X
rH = Y./X ;
end
end
end

7
Classes/01_transmit/PAMmapper.m
View File

@@ -24,7 +24,9 @@ classdef PAMmapper
if isa(signal_in,'Signal')
signal_in.signal = obj.map_(signal_in.signal);
signal_in = signal_in.normalize("mode","rms");
% signal_in = signal_in.normalize("mode","rms");
signal_in = signal_in.logbookentry();
out = signal_in;
else
@@ -42,7 +44,7 @@ classdef PAMmapper
function pam_sig = map_(obj,bitpattern)
switch obj.M
case 1
case 2
% 2-ASK: BPSK / OOK
pam_sig=bitpattern(:,1);
@@ -87,6 +89,7 @@ classdef PAMmapper
end
pam_sig = pam_sig/sqrt(21);
case 16
% 16-ASK:
x1 = bitpattern(:,1);

9
Classes/01_transmit/PAMsource.m
View File

@@ -9,6 +9,8 @@ classdef PAMsource
fsym
randkey
db_precode
mrds_code
mrds_blocklength
@@ -33,6 +35,8 @@ classdef PAMsource
options.fsym = 112e9;
options.randkey = 0;
options.db_precode = 0;
options.mrds_code = 0;
options.mrds_blocklength = 512;
@@ -90,6 +94,11 @@ classdef PAMsource
symbols = PAMmapper(obj.M,0).map(bits);
symbols.fs = obj.fsym;
if obj.db_precode
symbols = Duobinary().precode(symbols);
symbols = Duobinary().encode(symbols);
end
if obj.mrds_code
symbols = MRDS_coding("blocklength",obj.mrds_blocklength).encode(symbols);

23
Classes/04_DSP/A1_scheme.m
View File

@@ -1,23 +0,0 @@
classdef A1_scheme
%A1_SCHEME Summary of this class goes here
% Detailed explanation goes here
properties
Property1
end
methods
function obj = A1_scheme(inputArg1,inputArg2)
%A1_SCHEME Construct an instance of this class
% Detailed explanation goes here
obj.Property1 = inputArg1 + inputArg2;
end
function outputArg = method1(obj,inputArg)
%METHOD1 Summary of this method goes here
% Detailed explanation goes here
outputArg = obj.Property1 + inputArg;
end
end
end

245
Classes/04_DSP/Coding/Duobinary.m Normal file
View File

@@ -0,0 +1,245 @@
classdef Duobinary
%Duobinary Coding
% should work
properties(Access=public)
end
methods (Access=public)
function obj = Duobinary()
%NAME Construct an instance of this class
% Detailed explanation goes here
arguments
end
% %
% fn = fieldnames(options);
% for n = 1:numel(fn)
% try
% obj.(fn{n}) = options.(fn{n});
% end
% end
% do more stuff
end
function signalclass = precode(~,signalclass)
data = signalclass.signal;
u = unique(data);
M = numel(u);
%make unipolar
if M == 4
data = data .* sqrt(5);
elseif M == 6
data = data .* sqrt(10);
elseif M == 8
data = data .* sqrt(21);
elseif M == 16
data = data .* sqrt(85);
warning('Check if PAM16 implementation, mapping and scaling is correct!')
end
data = round(data);
b = min(data);
data = data - b;
data = data ./ 2;
assert(isequal((0:M-1)',unique(data)),'Check Duobinary Precoding'); %seems the signal is not unipolar
% Pre coding
bk = zeros(numel(data),1);
for k = 1:numel(data)-1
bk(k+1) =mod(data(k)-bk(k),M);
end
%make bipolar
bk = bk .* 2;
bk = bk + b;
if M == 4
bk = bk ./ sqrt(5);
elseif M == 6
bk = bk ./ sqrt(10);
elseif M == 8
bk = bk ./ sqrt(21);
end
signalclass.signal = bk;
assert(isequal(unique(signalclass.signal),u),'Check Duobinary Precoding'); %seems the signal is not the same as before
end
function signalclass = encode(~,signalclass)
data = signalclass.signal;
u = unique(data);
M = numel(u);
%make unipolar
if M == 4
data = data .* sqrt(5);
elseif M == 6
data = data .* sqrt(10);
elseif M == 8
data = data .* sqrt(21);
elseif M == 16
data = data .* sqrt(85);
warning('Check if PAM16 implementation, mapping and scaling is correct!')
end
data = round(data);
b = min(data);
data = data - b;
data = data ./ 2;
assert(isequal((0:M-1)',unique(data)),'Check Duobinary Precoding'); %seems the signal is not unipolar
% duobinary coding (1+D)
coeff = [1,1];
data = conv(data,coeff,"same");
%make bipolar
data = (data-round(mean(data),1));
if M == 4
data = data ./ sqrt(2.5); % 7-level constellation weighted with probability after DB code i.e. mean([-3 3 -2 -2 2 2 -1 -1 -1 1 1 1 0 0 0 0].^2) = 2.5 --> sqrt(2.5) == rms(constellation)
elseif M == 6
data = data ./ sqrt(5.8);
elseif M == 8
data = data ./ sqrt(10.5); % 15-level constellation weighted with probability after DB code i.e.
end
signalclass.signal = data;
% Code to evaluate the s in sqrt(s):
% M=6;
% c = (0:(2*M)-2) - ((2*M)-2)/2;
% s = 0;
% for i = 0:(2*M)-2
% p(i+1) = M-abs(i-(M-1));
% s = s+p(i+1)*c(i+1).^2;
% end
% s = s/M^2;
end
function signalclass = decode(~,signalclass)
data = signalclass.signal;
u = unique(data);
I = numel(u); %number of duobinary coded const. points
M = (I+1)/2; %PAM-M order
%make unipolar
if I == 7
data = data .* sqrt(2.5);
elseif I == 11
%todo
data = data .* sqrt(5.8);
warning('Check if PAM16 implementation, mapping and scaling is correct!')
elseif I == 15
data = data .* sqrt(10.5);
elseif I == 16
warning('Check if PAM16 implementation, mapping and scaling is correct!')
end
data = round(data);
b = min(data);
data = data - b;
data = round(data);
data = mod(data,M);
%make bipolar
data = data .* 2;
data = data - round(mean(data));
if M == 4
data = data ./ sqrt(5);
elseif M == 6
data = data ./ sqrt(10);
elseif M == 8
data = data ./ sqrt(21);
end
signalclass.signal = data;
end
function signalclass = feedbackdetector(~,signalclass)
data = signalclass.signal;
u = unique(data);
I = numel(u); %number of duobinary coded const. points
M = (I+1)/2; %PAM-M order
%make unipolar
if I == 7
data = data .* sqrt(2.5);
elseif I == 11
%todo
data = data .* sqrt(5.8);
warning('Check if PAM16 implementation, mapping and scaling is correct!')
elseif I == 15
data = data .* sqrt(10.5);
elseif I == 16
warning('Check if PAM16 implementation, mapping and scaling is correct!')
end
data = round(data);
b = min(data);
data = data - b;
data = round(data);
data_out = zeros(length(data),1);
const = 0:M-1;
%%FALSCH!
for k = 1:length(data)-1
d_ = data(k+1) - data_out(k);
[~,b] = min(abs(d_-[0:M-1]));
data_out(k+1) = const(b);
end
%make bipolar
data_out = data_out .* 2;
data_out = data_out - round(mean(data_out));
if M == 4
data_out = data_out ./ sqrt(5);
elseif M == 6
data_out = data_out ./ sqrt(10);
elseif M == 8
data_out = data_out ./ sqrt(21);
end
signalclass.signal = data_out;
end
end
methods (Access=private)
% Cant be seen from outside! So put all your functions here that can/
% shall not be called from outside
end
end

0
Classes/04_DSP/MRDS_coding.m → Classes/04_DSP/Coding/MRDS_coding.m
View File

0
Classes/04_DSP/EQ.m → Classes/04_DSP/Equalizer/EQ.m
View File

0
Classes/04_DSP/EQ_copy.m → Classes/04_DSP/Equalizer/EQ_copy.m
View File

972
Classes/04_DSP/EQ_silas.m → Classes/04_DSP/Equalizer/EQ_silas.m
View File

@@ -1,486 +1,486 @@
classdef EQ_silas < handle
%EQ_SILAS FFE and DFE Equalizer Playground
properties
% Important Signals
x_in %Input Sequence to be equalized
x_length
x_norm
d %reference signal
d_norm
d_constellation %constellation points of the reference
y_out %equalizer output signal
d_out %decision output
% FFE coefficients always named with "e"
Ne
Ce %memory length FFE
Ie1 %Indice Combination of 1nd order FFE
Ie2 %Indice Combination of 2nd order FFE
Ie3 %Indice Combination of 3nd order FFE
e %coefficients for FFE
% DFE coefficients always named with "b"
Nb
Cb %memory length DFE
Ib1 %Indice Combination of 1nd order DFE
Ib2 %Indice Combination of 2nd order DFE
Ib3 %Indice Combination of 3nd order DFE
b %coefficients for DFE
error
e_ffe
e_dfe
e_dc
% coefficients
mu_dc_train
mu_ffe_train
mu_dfe_train
mu_dc_dd
mu_ffe_dd
mu_dfe_dd
mu_combined_dd % [1st order FFE, 2nd order FFE, 3rd order FFE, all orders DFE]
delay
trainlength
sps
trainloops
ddloops
eq_parallelization_blocklength % block lengt of EQ (until now, only the dc subtraction is affected by this)
eq_updatelatency % time in symbols until the calculated updates reach the signal again (until now, only the dc subtraction is affected by this)
eq_avg_blocklength
end
methods
function obj = EQ_silas(options)
%EQ_SILAS Construct an instance of this class
arguments(Input)
options.Ne = [50 5 0] %Number of FFE coefficients (1st, 2nd and 3rd order)
options.Nb = [30 5 3] %Number of DFE coefficients (1st, 2nd and 3rd order)
options.trainloops = 2;
options.trainlength = 4096;
options.ddloops = 2;
options.delay = 0;
options.sps = 2;
options.mu_dc_train = 0.01;
options.mu_ffe_train = 0.005;
options.mu_dfe_train = 0.005;
options.mu_dc_dd = 0.01;
options.mu_ffe_dd = [0.0004 0.0005 0.0006];
options.mu_dfe_dd = 0.0005;
options.eq_parallelization_blocklength = 1;
options.eq_updatelatency = 1;
options.eq_avg_blocklength = 0;
end
fn = fieldnames(options);
for n = 1:numel(fn)
obj.(fn{n}) = options.(fn{n});
end
% Generate helpful vectors and initialize the filters with
% correct length:
obj.Ce = obj.calcVNLEMemoryLength(obj.Ne);
[obj.Ie2,obj.Ie3] = obj.calcIndiceVectors(obj.Ne);
obj.e = zeros(sum(obj.Ce),1);
obj.Cb = obj.calcVNLEMemoryLength(obj.Nb);
[obj.Ib2,obj.Ib3] = obj.calcIndiceVectors(obj.Nb);
obj.b = zeros(sum(obj.Cb),1);
end
function [signalclass_out,symbols_out] = process(obj,signalclass_in, reference_signalclass_in)
% actual processing of the signal (steps 1. - 3.)
% 1 normalize RMS
signalclass_in = signalclass_in.normalize("mode","rms");
% Process the EQ optimization
obj.process_(signalclass_in.signal', reference_signalclass_in.signal');
signalclass_in.signal = obj.y_out';
%change sampling frequency of outgoing signal
signalclass_in.fs = reference_signalclass_in.fs;
% append to logbook
lbdesc = ['EQ von Silas ist gelaufen '];
signalclass_in = signalclass_in.logbookentry(lbdesc);
symbols_out = signalclass_in;
symbols_out.signal = obj.d_out;
% write to output
signalclass_out = signalclass_in;
end
function process_(obj,x_in,d_in)
% 1) prepare signals
obj.e_dc = mean(x_in);
% 1.1) Input Signal
obj.x_in = [zeros(1,floor(obj.Ne(1)/2)) x_in zeros(1,obj.Ne(1))];
obj.x_length = length(x_in);
obj.x_norm = obj.calcPowerNormalization(x_in);
% 1.2 Reference Signal // Constellation
obj.d = [zeros(1,obj.Nb(1)-1) d_in zeros(1,obj.Nb(1))];
obj.d_constellation = unique(d_in);
obj.d_norm = obj.calcPowerNormalization(d_in);
% 1.3 Training
obj.trainingMode();
% 1.4 Decision Directed Mode
obj.decisionDirectedMode();
end
%% Adaptive Equalization Modes
function trainingMode(obj)
dc_block = ones(obj.eq_parallelization_blocklength,1);
for tloop = 1:obj.trainloops
m = 1+obj.delay;
dc_cnt = 0;
for n = obj.sps*obj.delay+1:obj.sps:obj.sps*obj.trainlength
m = m+1;
dc_cnt = dc_cnt+1;
%get Sigal input vectors with correct length for VNLE
x_in_block = obj.x_in(obj.Ne(1)+n+(obj.sps-1):-1:n+obj.sps).';
x_in_vnle_format = obj.calcVNLENonlinVecs(x_in_block,obj.Ie2,obj.Ie3,obj.Ne,obj.x_norm);
%get Reference input vectors with correct length for VNLE
d_block = obj.d(obj.Nb(1)-obj.delay+m-2:-1:m-obj.delay-1).';
d_vnle_format = obj.calcVNLENonlinVecs(d_block,obj.Ib2,obj.Ib3,obj.Nb,obj.d_norm);
obj.e_ffe = obj.e.' * x_in_vnle_format;
obj.e_dfe = obj.b.' * d_vnle_format;
% Calculate the Error
obj.error = obj.e_dc + obj.e_ffe - obj.e_dfe - obj.d(obj.Nb(1)-1+m-obj.delay);
if obj.mu_ffe_train ~= 0
%update FFE coefficients with LMS
obj.e = obj.e - obj.error*conj(x_in_vnle_format)*obj.mu_ffe_train;
else
%update FFE coefficients with NLMS
obj.e = obj.e - obj.error*x_in_vnle_format/(x_in_vnle_format.'*x_in_vnle_format);
end
%update DFE coefficients with LMS
obj.b = obj.b + obj.mu_dfe_train*obj.error*d_vnle_format;
%update DC error
dc_block(dc_cnt) = obj.error .* obj.mu_dc_train;
if dc_cnt == obj.eq_parallelization_blocklength
obj.e_dc = obj.e_dc - mean(dc_block(dc_cnt));
dc_cnt = 0;
end
end
end
end
function decisionDirectedMode(obj)
%start the dd mode with coefficients from training
coeff = [obj.e;obj.b];
obj.e_dc = ones(obj.eq_updatelatency,1).*obj.e_dc;
dc_block = ones(obj.eq_parallelization_blocklength,1);
for ddloop = 1:obj.ddloops
m = 0;
dc_cnt = 0;
mu_mat = diag([ones(1,obj.Ce(1))*obj.mu_ffe_dd(1)... %1st order ffe
ones(1,obj.Ce(2))*obj.mu_ffe_dd(2)... %2nd order ffe
ones(1,obj.Ce(3))*obj.mu_ffe_dd(3)... %3rd order ffe
ones(1,sum(obj.Cb))*obj.mu_dfe_dd]); %all order dfe
y = zeros(1,floor(obj.x_length/obj.sps));
d_feedback = zeros(obj.Cb(1),1);
d_vnle = obj.calcVNLENonlinVecs(d_feedback,obj.Ib2,obj.Ib3,obj.Nb,obj.d_norm);
d_hat = NaN(length(obj.d),numel(obj.d_constellation));
lvl_err_1 = NaN(length(obj.d),numel(obj.d_constellation));
lvl_err_2 = NaN(length(obj.d),numel(obj.d_constellation));
subtracted_error =NaN(length(obj.d),numel(obj.d_constellation));
y_1= NaN(length(obj.d),numel(obj.d_constellation));
y_2= NaN(length(obj.d),numel(obj.d_constellation));
lvl_err_mov = NaN(obj.eq_avg_blocklength,numel(obj.d_constellation));
m_reg = 0;
for k = 1:obj.sps:obj.x_length
dc_cnt = dc_cnt+1;
m=m+1;
%get Sigal input vectors with correct length for VNLE
x = obj.x_in(obj.Ne(1)+k-1:-1:k).';
%bring this signal to "special" VNLE format
x_vnle = obj.calcVNLENonlinVecs(x,obj.Ie2,obj.Ie3,obj.Ne,obj.x_norm);
%combine FFE with DFE to one vector (cursor between the two sequences)
x_d = [x_vnle;-d_vnle];
%Apply filter
y(m) = x_d.'* coeff;
%Decision 1
[~,symbol_idx] = min(abs(y(m) - obj.d_constellation)); % decision for closest constellation point
d_hat(m,symbol_idx) = obj.d_constellation(symbol_idx);
y_1(m,symbol_idx) = y(m); % after 1st iteration
%1st Error between FFE & DFE filtered signal and Decision
obj.error(m) = y(m) - d_hat(m,symbol_idx);
% lvl_err_1(m,symbol_idx) = y(m) - obj.d(m+1);
%
% %write current error to buffer
% lvl_err_mov(:,symbol_idx) = circshift(lvl_err_mov(:,symbol_idx),1);
% lvl_err_mov(1,symbol_idx) = obj.error(m);
%
% %Subtract a weighted error from y -> then Decision 2
% err = mean(lvl_err_mov(:,symbol_idx),'omitnan');
%
% y(m) = y(m)-(obj.mu_dc_dd(symbol_idx)*err);
%
% subtracted_error(m,symbol_idx) = obj.mu_dc_dd(symbol_idx)*mean(lvl_err_mov(:,symbol_idx),'omitnan');
%
% y_2(m,symbol_idx) = y(m);
%
% [~,symbol_idx] = min(abs(y(m) - obj.d_constellation)); % decision 2 for closest constellation point
%
% d_hat(m,symbol_idx) = obj.d_constellation(symbol_idx);
%
% obj.error(m) = y(m) - d_hat(m,symbol_idx);
%
% lvl_err_2(m,symbol_idx) = y(m) - obj.d(m+1);
%Update FFE and DFE coefficients
coeff = coeff - (mu_mat * (obj.error(m) * conj(x_d)));
% Append new decision to decision feedback
if obj.Nb(1) > 0
%shift up one index
d_feedback(2:end) = d_feedback(1:end-1);
%replace 1st index with current estimation
d_feedback(1) = d_hat(m,symbol_idx);
%build memorylike VNLE version
d_vnle = obj.calcVNLENonlinVecs(d_feedback,obj.Ib2,obj.Ib3,obj.Nb,obj.d_norm);
end
end
end
%%
obj.y_out = (circshift( y.' ,-(obj.delay))).';
obj.d_out = d_hat(1:2:end);
% evm1 = mean(lvl_err_1,'omitnan');
%
% evm2 = mean(lvl_err_2,'omitnan');
%
% figure(112)
% stem(evm1,'LineStyle','--','Marker','square','LineWidth',1);
% hold on;
% stem(evm2,'LineStyle',':','Marker','v','LineWidth',1);
% figure(14)
% scatter(1:length(lvl_err_1),subtracted_error,1,'.')
%
% lvl_err___ = lvl_err_true(~isnan(lvl_err_true));
% %lvl_err___ = lvl_err___-mean(lvl_err___);
% coeffs = arburg(lvl_err___,1000);
% fs_in = 92e9;
% [h,w] = freqz(1,coeffs,length(lvl_err___),"whole",fs_in);
% h = fftshift(h./max(abs(h)));
% freq_vec = linspace(-fs_in/2,fs_in/2,length(h));
% figure(111)
% hold on
% plot(freq_vec.*1e-9,20*log10(h),'DisplayName','burg');
%
% spectrum_plot(y,92e9);
%
% d = 2^nextpow2(length(y)/16);
%
% figure(1111)
% hold on
% pwelch(y,hamming(d),d/2,d,92e9,"centered","power");
end
%% Functions needed During Adaption
function x_in_vnle_format = calcVNLENonlinVecs(~,x_in_block,I_2,I_3,N_,norm_)
% These are the second and third order input signal products of the VNLE EQ
% ∑ h1 x_in(k-n1) + ∑∑ h2 x_in(k-n1)*x_in(k-n2) + ∑∑∑ h3 x_in(k-n1)*x_in(k-n2)*x_in(k-n3)
l1=length(x_in_block);
l2=length(I_2);
l3=length(I_3);
final_length = l1+l2+l3;
x_in_vnle_format = zeros(final_length,1);
idx = l1;
x_in_vnle_format(1:idx) = x_in_block;
if N_(2) > 0
delta_2 = round((N_(1)-N_(2)) / 2);
input_vec_se = x_in_block(delta_2:end) / norm_(2); %TODO normalization step
% Extract columns from I_2
col1 = input_vec_se(I_2(:,1));
col2 = input_vec_se(I_2(:,2));
x2 = col1 .* col2;
x_in_vnle_format(idx+1:idx+l2) = x2;
end
if N_(3) > 0
delta_3 = round((N_(1)-N_(3))/2);
input_vec_th = x_in_block(delta_3:end) / norm_(3);
% Extract columns from I_3
col1 = input_vec_th(I_3(:,1));
col2 = input_vec_th(I_3(:,2));
col3 = input_vec_th(I_3(:,3));
% Perform matrix multiplication
x3 = col1 .* col2 .* col3;
idx = idx+l2;
x_in_vnle_format(idx+1:idx+l3) = x3;
end
end
%% Functions needed for Preparation
function [C] = calcVNLEMemoryLength(~,N)
%calculates the memory length of VNLE
C = zeros(size(N));
for o = 1:numel(N)
switch o
case 1
C(o) = N(o);
case 2
C(o) = N(o)*(N(o)+1) / 2;
case 3
C(o) = N(o)*(N(o)+1)*(N(o)+2) / 6;
end
end
end
function [indvec2nd, indvec3rd] = calcIndiceVectors(~,N)
% Init vectors of 2nd and 3rd order coefficient indices ->
% yield combination with
for order = 2:numel(N)
n = N(order);
v = 1:n; % Ursprünglicher Vektor
row = 1;
% Schleifen zur Generierung des Indize Vektors
switch order
case 2
indvec2nd = zeros(n*(n+1)/2, order);
for i = 1:n
for j = i:n
indvec2nd(row, :) = [v(i) v(j)];
row = row + 1;
end
end
case 3
indvec3rd = zeros(n*(n+1)*(n+2)/6, 3);
for i = 1:n
for j = i:n
for k = j:n
indvec3rd(row, :) = [v(i) v(j) v(k)];
row = row + 1;
end
end
end
end
end
end
function powerNorm = calcPowerNormalization(~,v)
powerNorm(1) = sqrt(mean(abs(v ).^2));
powerNorm(2) = sqrt(mean(abs(v.^2).^2));
powerNorm(3) = sqrt(mean(abs(v.^3).^2));
end
end
end
classdef EQ_silas < handle
%EQ_SILAS FFE and DFE Equalizer Playground
properties
% Important Signals
x_in %Input Sequence to be equalized
x_length
x_norm
d %reference signal
d_norm
d_constellation %constellation points of the reference
y_out %equalizer output signal
d_out %decision output
% FFE coefficients always named with "e"
Ne
Ce %memory length FFE
Ie1 %Indice Combination of 1nd order FFE
Ie2 %Indice Combination of 2nd order FFE
Ie3 %Indice Combination of 3nd order FFE
e %coefficients for FFE
% DFE coefficients always named with "b"
Nb
Cb %memory length DFE
Ib1 %Indice Combination of 1nd order DFE
Ib2 %Indice Combination of 2nd order DFE
Ib3 %Indice Combination of 3nd order DFE
b %coefficients for DFE
error
e_ffe
e_dfe
e_dc
% coefficients
mu_dc_train
mu_ffe_train
mu_dfe_train
mu_dc_dd
mu_ffe_dd
mu_dfe_dd
mu_combined_dd % [1st order FFE, 2nd order FFE, 3rd order FFE, all orders DFE]
delay
trainlength
sps
trainloops
ddloops
eq_parallelization_blocklength % block lengt of EQ (until now, only the dc subtraction is affected by this)
eq_updatelatency % time in symbols until the calculated updates reach the signal again (until now, only the dc subtraction is affected by this)
eq_avg_blocklength
end
methods
function obj = EQ_silas(options)
%EQ_SILAS Construct an instance of this class
arguments(Input)
options.Ne = [50 5 0] %Number of FFE coefficients (1st, 2nd and 3rd order)
options.Nb = [30 5 3] %Number of DFE coefficients (1st, 2nd and 3rd order)
options.trainloops = 2;
options.trainlength = 4096;
options.ddloops = 2;
options.delay = 0;
options.sps = 2;
options.mu_dc_train = 0.01;
options.mu_ffe_train = 0.005;
options.mu_dfe_train = 0.005;
options.mu_dc_dd = 0.01;
options.mu_ffe_dd = [0.0004 0.0005 0.0006];
options.mu_dfe_dd = 0.0005;
options.eq_parallelization_blocklength = 1;
options.eq_updatelatency = 1;
options.eq_avg_blocklength = 0;
end
fn = fieldnames(options);
for n = 1:numel(fn)
obj.(fn{n}) = options.(fn{n});
end
% Generate helpful vectors and initialize the filters with
% correct length:
obj.Ce = obj.calcVNLEMemoryLength(obj.Ne);
[obj.Ie2,obj.Ie3] = obj.calcIndiceVectors(obj.Ne);
obj.e = zeros(sum(obj.Ce),1);
obj.Cb = obj.calcVNLEMemoryLength(obj.Nb);
[obj.Ib2,obj.Ib3] = obj.calcIndiceVectors(obj.Nb);
obj.b = zeros(sum(obj.Cb),1);
end
function [signalclass_out,symbols_out] = process(obj,signalclass_in, reference_signalclass_in)
% actual processing of the signal (steps 1. - 3.)
% 1 normalize RMS
signalclass_in = signalclass_in.normalize("mode","rms");
% Process the EQ optimization
obj.process_(signalclass_in.signal', reference_signalclass_in.signal');
signalclass_in.signal = obj.y_out';
%change sampling frequency of outgoing signal
signalclass_in.fs = reference_signalclass_in.fs;
% append to logbook
lbdesc = ['EQ von Silas ist gelaufen '];
signalclass_in = signalclass_in.logbookentry(lbdesc);
symbols_out = signalclass_in;
symbols_out.signal = obj.d_out;
% write to output
signalclass_out = signalclass_in;
end
function process_(obj,x_in,d_in)
% 1) prepare signals
obj.e_dc = mean(x_in);
% 1.1) Input Signal
obj.x_in = [zeros(1,floor(obj.Ne(1)/2)) x_in zeros(1,obj.Ne(1))];
obj.x_length = length(x_in);
obj.x_norm = obj.calcPowerNormalization(x_in);
% 1.2 Reference Signal // Constellation
obj.d = [zeros(1,obj.Nb(1)-1) d_in zeros(1,obj.Nb(1))];
obj.d_constellation = unique(d_in);
obj.d_norm = obj.calcPowerNormalization(d_in);
% 1.3 Training
obj.trainingMode();
% 1.4 Decision Directed Mode
obj.decisionDirectedMode();
end
%% Adaptive Equalization Modes
function trainingMode(obj)
dc_block = ones(obj.eq_parallelization_blocklength,1);
for tloop = 1:obj.trainloops
m = 1+obj.delay;
dc_cnt = 0;
for n = obj.sps*obj.delay+1:obj.sps:obj.sps*obj.trainlength
m = m+1;
dc_cnt = dc_cnt+1;
%get Sigal input vectors with correct length for VNLE
x_in_block = obj.x_in(obj.Ne(1)+n+(obj.sps-1):-1:n+obj.sps).';
x_in_vnle_format = obj.calcVNLENonlinVecs(x_in_block,obj.Ie2,obj.Ie3,obj.Ne,obj.x_norm);
%get Reference input vectors with correct length for VNLE
d_block = obj.d(obj.Nb(1)-obj.delay+m-2:-1:m-obj.delay-1).';
d_vnle_format = obj.calcVNLENonlinVecs(d_block,obj.Ib2,obj.Ib3,obj.Nb,obj.d_norm);
obj.e_ffe = obj.e.' * x_in_vnle_format;
obj.e_dfe = obj.b.' * d_vnle_format;
% Calculate the Error
obj.error = obj.e_dc + obj.e_ffe - obj.e_dfe - obj.d(obj.Nb(1)-1+m-obj.delay);
if obj.mu_ffe_train ~= 0
%update FFE coefficients with LMS
obj.e = obj.e - obj.error*conj(x_in_vnle_format)*obj.mu_ffe_train;
else
%update FFE coefficients with NLMS
obj.e = obj.e - obj.error*x_in_vnle_format/(x_in_vnle_format.'*x_in_vnle_format);
end
%update DFE coefficients with LMS
obj.b = obj.b + obj.mu_dfe_train*obj.error*d_vnle_format;
%update DC error
dc_block(dc_cnt) = obj.error .* obj.mu_dc_train;
if dc_cnt == obj.eq_parallelization_blocklength
obj.e_dc = obj.e_dc - mean(dc_block(dc_cnt));
dc_cnt = 0;
end
end
end
end
function decisionDirectedMode(obj)
%start the dd mode with coefficients from training
coeff = [obj.e;obj.b];
obj.e_dc = ones(obj.eq_updatelatency,1).*obj.e_dc;
dc_block = ones(obj.eq_parallelization_blocklength,1);
for ddloop = 1:obj.ddloops
m = 0;
dc_cnt = 0;
mu_mat = diag([ones(1,obj.Ce(1))*obj.mu_ffe_dd(1)... %1st order ffe
ones(1,obj.Ce(2))*obj.mu_ffe_dd(2)... %2nd order ffe
ones(1,obj.Ce(3))*obj.mu_ffe_dd(3)... %3rd order ffe
ones(1,sum(obj.Cb))*obj.mu_dfe_dd]); %all order dfe
y = zeros(1,floor(obj.x_length/obj.sps));
d_feedback = zeros(obj.Cb(1),1);
d_vnle = obj.calcVNLENonlinVecs(d_feedback,obj.Ib2,obj.Ib3,obj.Nb,obj.d_norm);
d_hat = NaN(length(obj.d),numel(obj.d_constellation));
lvl_err_1 = NaN(length(obj.d),numel(obj.d_constellation));
lvl_err_2 = NaN(length(obj.d),numel(obj.d_constellation));
subtracted_error =NaN(length(obj.d),numel(obj.d_constellation));
y_1= NaN(length(obj.d),numel(obj.d_constellation));
y_2= NaN(length(obj.d),numel(obj.d_constellation));
lvl_err_mov = NaN(obj.eq_avg_blocklength,numel(obj.d_constellation));
m_reg = 0;
for k = 1:obj.sps:obj.x_length
dc_cnt = dc_cnt+1;
m=m+1;
%get Sigal input vectors with correct length for VNLE
x = obj.x_in(obj.Ne(1)+k-1:-1:k).';
%bring this signal to "special" VNLE format
x_vnle = obj.calcVNLENonlinVecs(x,obj.Ie2,obj.Ie3,obj.Ne,obj.x_norm);
%combine FFE with DFE to one vector (cursor between the two sequences)
x_d = [x_vnle;-d_vnle];
%Apply filter
y(m) = x_d.'* coeff;
%Decision 1
[~,symbol_idx] = min(abs(y(m) - obj.d_constellation)); % decision for closest constellation point
d_hat(m,symbol_idx) = obj.d_constellation(symbol_idx);
y_1(m,symbol_idx) = y(m); % after 1st iteration
%1st Error between FFE & DFE filtered signal and Decision
obj.error(m) = y(m) - d_hat(m,symbol_idx);
% lvl_err_1(m,symbol_idx) = y(m) - obj.d(m+1);
%
% %write current error to buffer
% lvl_err_mov(:,symbol_idx) = circshift(lvl_err_mov(:,symbol_idx),1);
% lvl_err_mov(1,symbol_idx) = obj.error(m);
%
% %Subtract a weighted error from y -> then Decision 2
% err = mean(lvl_err_mov(:,symbol_idx),'omitnan');
%
% y(m) = y(m)-(obj.mu_dc_dd(symbol_idx)*err);
%
% subtracted_error(m,symbol_idx) = obj.mu_dc_dd(symbol_idx)*mean(lvl_err_mov(:,symbol_idx),'omitnan');
%
% y_2(m,symbol_idx) = y(m);
%
% [~,symbol_idx] = min(abs(y(m) - obj.d_constellation)); % decision 2 for closest constellation point
%
% d_hat(m,symbol_idx) = obj.d_constellation(symbol_idx);
%
% obj.error(m) = y(m) - d_hat(m,symbol_idx);
%
% lvl_err_2(m,symbol_idx) = y(m) - obj.d(m+1);
%Update FFE and DFE coefficients
coeff = coeff - (mu_mat * (obj.error(m) * conj(x_d)));
% Append new decision to decision feedback
if obj.Nb(1) > 0
%shift up one index
d_feedback(2:end) = d_feedback(1:end-1);
%replace 1st index with current estimation
d_feedback(1) = d_hat(m,symbol_idx);
%build memorylike VNLE version
d_vnle = obj.calcVNLENonlinVecs(d_feedback,obj.Ib2,obj.Ib3,obj.Nb,obj.d_norm);
end
end
end
%%
obj.y_out = (circshift( y.' ,-(obj.delay))).';
obj.d_out = d_hat(1:2:end);
% evm1 = mean(lvl_err_1,'omitnan');
%
% evm2 = mean(lvl_err_2,'omitnan');
%
% figure(112)
% stem(evm1,'LineStyle','--','Marker','square','LineWidth',1);
% hold on;
% stem(evm2,'LineStyle',':','Marker','v','LineWidth',1);
% figure(14)
% scatter(1:length(lvl_err_1),subtracted_error,1,'.')
%
% lvl_err___ = lvl_err_true(~isnan(lvl_err_true));
% %lvl_err___ = lvl_err___-mean(lvl_err___);
% coeffs = arburg(lvl_err___,1000);
% fs_in = 92e9;
% [h,w] = freqz(1,coeffs,length(lvl_err___),"whole",fs_in);
% h = fftshift(h./max(abs(h)));
% freq_vec = linspace(-fs_in/2,fs_in/2,length(h));
% figure(111)
% hold on
% plot(freq_vec.*1e-9,20*log10(h),'DisplayName','burg');
%
% spectrum_plot(y,92e9);
%
% d = 2^nextpow2(length(y)/16);
%
% figure(1111)
% hold on
% pwelch(y,hamming(d),d/2,d,92e9,"centered","power");
end
%% Functions needed During Adaption
function x_in_vnle_format = calcVNLENonlinVecs(~,x_in_block,I_2,I_3,N_,norm_)
% These are the second and third order input signal products of the VNLE EQ
% ∑ h1 x_in(k-n1) + ∑∑ h2 x_in(k-n1)*x_in(k-n2) + ∑∑∑ h3 x_in(k-n1)*x_in(k-n2)*x_in(k-n3)
l1=length(x_in_block);
l2=length(I_2);
l3=length(I_3);
final_length = l1+l2+l3;
x_in_vnle_format = zeros(final_length,1);
idx = l1;
x_in_vnle_format(1:idx) = x_in_block;
if N_(2) > 0
delta_2 = round((N_(1)-N_(2)) / 2);
input_vec_se = x_in_block(delta_2:end) / norm_(2); %TODO normalization step
% Extract columns from I_2
col1 = input_vec_se(I_2(:,1));
col2 = input_vec_se(I_2(:,2));
x2 = col1 .* col2;
x_in_vnle_format(idx+1:idx+l2) = x2;
end
if N_(3) > 0
delta_3 = round((N_(1)-N_(3))/2);
input_vec_th = x_in_block(delta_3:end) / norm_(3);
% Extract columns from I_3
col1 = input_vec_th(I_3(:,1));
col2 = input_vec_th(I_3(:,2));
col3 = input_vec_th(I_3(:,3));
% Perform matrix multiplication
x3 = col1 .* col2 .* col3;
idx = idx+l2;
x_in_vnle_format(idx+1:idx+l3) = x3;
end
end
%% Functions needed for Preparation
function [C] = calcVNLEMemoryLength(~,N)
%calculates the memory length of VNLE
C = zeros(size(N));
for o = 1:numel(N)
switch o
case 1
C(o) = N(o);
case 2
C(o) = N(o)*(N(o)+1) / 2;
case 3
C(o) = N(o)*(N(o)+1)*(N(o)+2) / 6;
end
end
end
function [indvec2nd, indvec3rd] = calcIndiceVectors(~,N)
% Init vectors of 2nd and 3rd order coefficient indices ->
% yield combination with
for order = 2:numel(N)
n = N(order);
v = 1:n; % Ursprünglicher Vektor
row = 1;
% Schleifen zur Generierung des Indize Vektors
switch order
case 2
indvec2nd = zeros(n*(n+1)/2, order);
for i = 1:n
for j = i:n
indvec2nd(row, :) = [v(i) v(j)];
row = row + 1;
end
end
case 3
indvec3rd = zeros(n*(n+1)*(n+2)/6, 3);
for i = 1:n
for j = i:n
for k = j:n
indvec3rd(row, :) = [v(i) v(j) v(k)];
row = row + 1;
end
end
end
end
end
end
function powerNorm = calcPowerNormalization(~,v)
powerNorm(1) = sqrt(mean(abs(v ).^2));
powerNorm(2) = sqrt(mean(abs(v.^2).^2));
powerNorm(3) = sqrt(mean(abs(v.^3).^2));
end
end
end

0
Classes/04_DSP/EQ_silas_ofc.m → Classes/04_DSP/Equalizer/EQ_silas_ofc.m
View File

0
Classes/04_DSP/EQ_silas_plain.m → Classes/04_DSP/Equalizer/EQ_silas_plain.m
View File

0
Classes/04_DSP/EQ_silas_sliding_window_dc_removal.m → Classes/04_DSP/Equalizer/EQ_silas_sliding_window_dc_removal.m
View File

0
Classes/04_DSP/FFE.m → Classes/04_DSP/Equalizer/FFE.m
View File

0
Classes/04_DSP/FFE_DCremoval.m → Classes/04_DSP/Equalizer/FFE_DCremoval.m
View File

0
Classes/04_DSP/FFE_DFE.m → Classes/04_DSP/Equalizer/FFE_DFE.m
View File

0
Classes/04_DSP/FFE_FFDCAVG.m → Classes/04_DSP/Equalizer/FFE_FFDCAVG.m
View File

0
Classes/04_DSP/FFE_adaptive_decision.m → Classes/04_DSP/Equalizer/FFE_adaptive_decision.m
View File

0
Classes/04_DSP/VNLE.m → Classes/04_DSP/Equalizer/VNLE.m
View File

8
Datatypes/delay_mode.m Normal file
View File

@@ -0,0 +1,8 @@
classdef delay_mode < int32
enumeration
samples (1)
time (2)
end
end

4
Functions/calc_ber.m
View File

@@ -10,6 +10,8 @@ end
options.skip_end = abs(options.skip_end);
options.skip_front = abs(options.skip_front);
assert((options.skip_end+options.skip_front)<length(data_in),"You can not skip more bits than overall length of data! Set skip_front or skip_end to lower value or check data_in");
bits = 0;
errors= 0;
ber= 0;
@@ -29,7 +31,7 @@ if length(data_ref) == length(data_in)
else
errorIndice = sum(data_in ~= data_ref,1);
errors = sum(errorIndice ,"all" );
[~,errorIndice] = find(errorIndice==1);
[~,errorIndice] = find(errorIndice~=0);
end
catch

BIN
projects/standard_system.mat Normal file
View File

Binary file not shown.

BIN
projects/standard_system/400G_simulative_setup.mat Normal file
View File

Binary file not shown.

BIN
projects/standard_system/bla.mat Normal file
View File

Binary file not shown.

28
projects/standard_system/freqresp_test.m Normal file
View File

@@ -0,0 +1,28 @@
measure = 0;
if measure
freqresp = ChannelFreqResp("Nacq",1024,"Navg",64,"Ncp",70,"f_ref",256e9);
%
Digi_sig = freqresp.buildOFDM();
else
[Digi_sig,Symbols,Bits] = PAMsource("fsym",fsym,"M",M,"order",18,"useprbs",0,...
"fs_out",M8199.fdac,"applyclipping",1,"clipfactor",1.7,"applypulseform",1,"pulseformer",Pform,"randkey",pn_key,"mrds_code",usemrds,"mrds_blocklength",512).process();
end
Digi_sig.spectrum("fignum",1112,"displayname",['Signal']);
maxamp = -1;
El_sig = freqresp.precomp(Digi_sig,"maxampdb",maxamp);
El_sig.spectrum("fignum",1112,"displayname",['maxamp:',num2str(maxamp)]);
El_sig = Filter('filtdegree',2,"f_cutoff",60e9,"fs",256e9,"filterType",filtertypes.butterworth,"active",true).process(El_sig);
if measure
freqresp.estimate(El_sig,"fileName",'','save',false);
end
El_sig.spectrum("fignum",1112,"displayname",['after filter; maxamp:',num2str(maxamp)]);

197
projects/standard_system/imdd_minimal_2.m Normal file
View File

@@ -0,0 +1,197 @@
%% Parameter to simulate and save
params = struct;
params.M = [4];
params.datarate = [440];
params.rop = [0];
precomp_mode = 2; %0=do nothing ; 1= measure; 2=precomp active
precomp_path = "C:\Users\Silas\Documents\MATLAB\imdd_simulation\projects\standard_system\";
precomp_fn = "400G_simulative_setup";
usemrds = 0;
name = ['wh_',strrep(num2str(now),'.','')];
wh = DataStorage(params);
wh.addStorage("ber_ffe");
%% Init Params
link_length = 0; %meter
pn_key = 2;
laser_linewidth = 0;
endcnt = prod(wh.dim);
cnt=0;
disp(['Start Simulation of ',num2str(endcnt),' loops...'])
tic
for M = wh.parameter.M.values
for datarate = wh.parameter.datarate.values
% SETUP HERE: %%
kover = 16;
M8199 = M8199B("kover",kover);
fdac = M8199.fdac;
fsym = round(datarate / log2(M)) * 1e9;
rrcalpha = 0.05;
Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"rrcalpha",rrcalpha);
% MAIN SIGNAL
%%%%% Symbol Generation %%%%%%
[Digi_sig,Symbols,Bits] = PAMsource("fsym",fsym,"M",M,"order",18,"useprbs",0,...
"fs_out",M8199.fdac,"applyclipping",1,"clipfactor",1.7,...
"applypulseform",0,"pulseformer",Pform,"randkey",pn_key,...
"db_precode",1,...
"mrds_code",usemrds,"mrds_blocklength",512).process();
Digi_sig.eye(fsym,M);
Digi_sig.spectrum("fignum",11,"displayname",'before precomp');
if precomp_mode == 1
freqresp = ChannelFreqResp("Nacq",1024,"Navg",64,"Ncp",63,'f_ref',Digi_sig.fs);
Digi_sig = freqresp.buildOFDM();
elseif precomp_mode == 2
Digi_sig = freqresp.precomp(Digi_sig,'maxampdb',3,'loadPath',precomp_path,'fileName',precomp_fn);
Digi_sig.spectrum("fignum",11,"displayname",'after precomp');
end
%%%%% AWG %%%%%%
El_sig = M8199.process(Digi_sig);
% El_sig.signal = awgn(El_sig.signal,-3,'measured',pn_key);
%%%%% Lowpass el. components %%%%%%
El_sig = Filter('filtdegree',2,"f_cutoff",60e9,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(El_sig);
%%%%% Electrical Driver Amplifier %%%%%%
El_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","gain","amplification_db",9).process(El_sig);
% El_sig = El_sig.setPower(6,"dBm");
fprintf('Driver output power: %s dBm\n', num2str(El_sig.power));
fprintf('Driver output peak voltage: %s Vpp \n', num2str(max(El_sig.signal)-min(El_sig.signal)));
% MAIN SIGNAL
%%%%% MODULATE E/O CONVERSION %%%%%%
vbias_rel = 0.5;
u_pi = 2.9;
vbias = -vbias_rel*u_pi;
[Opt_sig] = EML("mode",eml_mode.im_cosinus,"power",3,"fsimu",El_sig.fs,"lambda",1290,"bias",vbias,"u_pi",u_pi,"linewidth",laser_linewidth,"randomkey",pn_key).process(El_sig);
% Opt_sig.eye(fsym,M);
% %
% figure(10)
% hold on
% scatter(El_sig.signal(1:100000)+vbias,(abs(Opt_sig.signal(1:100000)).^2)*1e3,0.1,'.','DisplayName','Modulator TF')
% ylim([0 4]);
% xlim([-u_pi/2, u_pi/2]+vbias);
% xlabel('Input in V')
% ylabel('abs(Output) in mW')
Optfilter = Filter('filtdegree',6,"f_cutoff",fsym.*0.7,"fs",fdac*kover,"filterType",filtertypes.gaussian,"active",true);
Opt_sig = Optfilter.process(Opt_sig);
Opt_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",0).process(Opt_sig);
i_ = wh.parameter.rop.length;
ber_ffe=zeros(i_);
patten=zeros(i_);
%%%%% Interference Signal Fiber Prop %%%%%%
Opt_sig = Fiber("fsimu",Opt_sig.fs,"fiber_length",link_length/1000,"alpha",0.3,"D",0,"lambda0",1310,"gamma",0,"Dslope",0.07).process(Opt_sig);
% Receiver ROP curve
for i = 1:i_
rop=wh.parameter.rop.values(i);
% Set ROP
Rx_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",rop).process(Opt_sig);
patten(i) = Rx_sig.power;
%%%%%% Square Law %%%%%%
Rx_sig = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20,"nep",1.8e-11).process(Rx_sig);
%%%%%% Lowpass PhDiode %%%%%%
Rx_sig = Filter('filtdegree',2,"f_cutoff",70e9,"fs",fdac*kover,"filterType",filtertypes.gaussian,"active",true).process(Rx_sig);
%%%%%% Scope %%%%%%
fadc = 256e9;
Lp_scpe = Filter('filtdegree',4,"f_cutoff",100e9,"fs",fadc,"filterType",filtertypes.butterworth,"active",true);
Scpe_sig = Scope("fsimu",fdac*kover,"fadc",fadc,...
"delay",0,"fixed_delay",0,"filtertype",filtertypes.butterworth,...
"samplingdelay",0,"rand_samplingdelay",0,"freq_offset",0,"samp_jitter",0,...
"adcresolution",10,"quantbuffer",0.1,'block_dc',1,'lpf_active',1,'H_lpf',Lp_scpe).process(Rx_sig);
if precomp_mode == 1
freqresp.estimate(Scpe_sig,"save",true,"savePath",precomp_path,"fileName",precomp_fn);
freqresp.plot();
end
%%%%%% Sample to 2x fsym %%%%%%
Scpe_sig = Scpe_sig.resample("fs_in",fadc,"fs_out",2*fsym);
%%%%%% Sync Rx signal with reference %%%%%%
[Scpe_sig,D,cuts] = Scpe_sig.tsynch("reference",Symbols,"fs_ref",fsym);
%%%%% EQUALIZE %%%%%%
%
Eq = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",25,"sps",2,"decide",1);
%Eq = VNLE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",[0.0004 0.0005 0.0006],"mu_tr",0,"order",[50,7,7],"sps",2,"decide",1);
[EQ_sig] = Eq.process(Scpe_sig,Symbols);
Rx_bits = PAMmapper(M,0).demap(EQ_sig);
[~,errors_bm,ber_ffe(i),errors] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
disp(['BER: ',sprintf('%.1E',ber_ffe(i)),' - - ROP: ',num2str(patten(i)),'dBm - - PAM-',num2str(M),' - - ',num2str(fsym*1e-9),' GBd']);
end
for i = 1:i_
rop=wh.parameter.rop.values(i);
wh.addValueToStorage(ber_ffe(i) ,'ber_ffe',M,datarate,rop);
end
toc
disp(['Simulated: ',num2str(cnt/endcnt*100),' %']);
wh.save('C:\Users\Silas\Documents\MATLAB\imdd_simulation\projects\MPI_August\auswertung\')
end
end
cols = linspecer(8);
%cnt = cnt+1;
ber_ffe = wh.getStoValue('ber_ffe',M,datarate,wh.parameter.rop.values);
% Create the initial plot
figure(44);
a = gca;
hold on; % Retain the plot so new points can be added without complete redraw
plot(wh.parameter.rop.values,ber_ffe,"LineWidth",0.5,"LineStyle","-","Marker",".","MarkerSize",15,"DisplayName","FFE only");
yline(3.8e-3,'DisplayName','HD-FEC','LineStyle','--','HandleVisibility','off');
xlabel('Received Optical Power (dBm)');
ylabel('Bit Error Rate (BER)');
title('Bit Error Rate vs. ROP');
set(gca,'yscale','log');
set(gca,'Box','on');
grid on;
grid minor
legend

54
projects/standard_system/sensitivity_analysis.m Normal file
View File

@@ -0,0 +1,54 @@
M = wh.parameter.M.values(1);
datarates = wh.parameter.datarate.values;
rop = wh.parameter.rop.values;
cnt = 1;
for dr = datarates
ber_ffe = wh.getStoValue('ber_ffe',M,dr,rop);
sens_ffe(cnt) = getIntersection(ber_ffe,rop);
figure(112)
hold on
plot(rop,ber_ffe,"LineWidth",0.5,"LineStyle","-","Marker",".","MarkerSize",15,"DisplayName",['FFE only; ',num2str(dr), ' Gbps']);
yline(3.8e-3,'DisplayName','HD-FEC','LineStyle','--','HandleVisibility','off');
xlabel('Received Optical Power (dBm)');
ylabel('Bit Error Rate (BER)');
title('Bit Error Rate vs. ROP');
set(gca,'yscale','log');
set(gca,'Box','on');
grid on;
grid minor
legend
cnt = cnt+1;
end
figure(223);
hold on; % Retain the plot so new points can be added without complete redraw
plot(datarates,sens_ffe,"LineWidth",1,"LineStyle","-","Marker",".","MarkerSize",10,"DisplayName",['Linewidth: ',num2str(laser_linewidth.*1e-6),' MHz']);
xlabel('Signal to Interference Ratio (dB)');
ylabel('Receiver Sensitivity');
title(['Bit Error Rate vs. SIR']);
grid on;
legend
function i = getIntersection(ber,rop)
%get intersection between rop curve and hd-fec limit
hdfec = 3.8e-3 .* ones(size(ber));
i = InterX([rop;hdfec'],[rop;ber']);
if isempty(i)
i = NaN;
else
i = i(1);
end
end

36
projects/standard_system/testcode_duobinary.m Normal file
View File

@@ -0,0 +1,36 @@
% Define parameters
n = 1; % Define the memory length n+1 (e.g., 1 for duobinary, 2 for tribinary, 3 for tetrabinary)
M = 4; % Modulation order (e.g., 4 for QPSK)
d = randi([0 M-1], 1, 100); % Example input sequence d(k)
% Pre-coding: bPR(k) = (d(k) - bPR(k-1)) mod M
bPR = zeros(size(d));
%bPR(1) = d(1); % Initial condition for pre-coding
for k = 1:numel(d)-1
bPR(k+1) = mod( d(k) - bPR(k), M );
end
% Polybinary coding: Convolution with (1 + z^-1)^n
% Generate coefficients for the polybinary filter (1 + z^-1)^n
coeff = ones(1, n+1); % Initialize coefficients with ones
for i = 2:n+1
coeff(i) = nchoosek(n, i-1); % Binomial coefficient for expansion
end
% Apply convolution to get the polybinary coded sequence
polybinary_data = conv(bPR, coeff, "same");
% Modulo operation for decoding to stay within M levels
d_reconstructed = mod(polybinary_data, M);
% Check if reconstruction is correct
reconstruction_success = isequal(d, d_reconstructed);
figure(111)
hold on
stairs(d_reconstructed)
stairs(d)
disp(['Reconstruction success: ', num2str(reconstruction_success)]);

BIN
projects/standard_system/tf_test.mat Normal file
View File

Binary file not shown.

55
test/duobinary_minimal_example.m Normal file
View File

@@ -0,0 +1,55 @@
useprbs = 1;
M = 8;
randkey = 1;
%%%%% PRBS Generation in correct shape for Modulation Format %%%%%%
O = 18; %order of prbs
N = 2^(O-1); %length of prbs
[~,seed] = prbs(O,1); %initialize first seed of prbs
bitpattern=[];
if useprbs
for i = 1:log2(M)
[bitpattern(:,i),seed] = prbs(O,N,seed);
end
else
s = RandStream('twister','Seed',randkey);
for i = 1:log2(M)
bitpattern(:,i) = randi(s,[0 1], N, 1);
end
end
if M == 6
bitpattern = reshape(bitpattern,[],1);
bitpattern = bitpattern(1:end-mod(length(bitpattern),5));
end
Tx_bits = Informationsignal(bitpattern);
Symbols_tx = PAMmapper(M,0).map(Tx_bits);
Symbols_tx.fs = fsym;
% Symbols = Duobinary().precode(Symbols_tx);
Symbols = Duobinary().encode(Symbols_tx);
% Symbols.spectrum("displayname","db","fignum",12)
Symbols = Duobinary().feedbackdetector(Symbols);
% Symbols.eye(fsym,M);
% Symbols = Duobinary().decode(Symbols);
Rx_bits = PAMmapper(M,0).demap(Symbols);
[~,error_num,ber,error_pos] = calc_ber(Tx_bits.signal,Rx_bits.signal,"skip_front",0,"skip_end",1,"returnErrorLocation",1);
disp(['BER: ',sprintf('%.1E',ber),' - - PAM-',num2str(M)]);
%
figure;hold on
stairs(Rx_bits.signal(1:100,1),'LineWidth',2)
stairs(Tx_bits.signal(1:100,1),'LineStyle',':','LineWidth',2);
figure;hold on
stairs(Symbols_tx.signal(1:100,1),'LineWidth',2)
stairs(Symbols.signal(1:100,1),'LineStyle',':','LineWidth',2);