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

MPI Simulations and stuff

Browse Source
This commit is contained in:
Silas Oettinghaus
2024-08-14 09:36:51 +02:00
parent 1eeb970d8f
commit 34f9149346
61 changed files with 4295 additions and 429 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
Classes/00_signals/Electricalsignal.m
View File

@@ -3,7 +3,7 @@ classdef Electricalsignal < Signal
% Detailed explanation goes here
properties
fs
end
methods
@@ -26,7 +26,7 @@ classdef Electricalsignal < Signal
end
function pow = power(obj)
function pow = power50(obj)
% Power of an electrical signal
R = 50;

BIN
Classes/00_signals/Frame 1.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 101 KiB

12
Classes/00_signals/Informationsignal.m
View File

@@ -3,7 +3,7 @@ classdef Informationsignal < Signal
% Detailed explanation goes here
properties
fs
end
methods
@@ -26,11 +26,11 @@ classdef Informationsignal < Signal
end
function pow = power(obj)
pow = mean(abs(obj.signal.^2),"all") ;
end
% function pow = power(obj)
%
% pow = mean(abs(obj.signal.^2),"all") ;
%
% end
end

19
Classes/00_signals/Opticalsignal.m
View File

@@ -5,7 +5,6 @@ classdef Opticalsignal < Signal
properties
nase
lambda
fs
end
@@ -52,14 +51,6 @@ classdef Opticalsignal < Signal
end
function pow = power(obj)
pow = mean( abs(obj.signal).^2 ) ;
pow = pow2db(pow)+30; %dbm
end
function cspr = cspr(obj)
carrier_power_dbm = pow2db( abs(mean(obj.signal)).^2 )+30; % dB -> +30 -> dBm
@@ -70,17 +61,15 @@ classdef Opticalsignal < Signal
cspr_lin = pow2db(carr_lin/sign_lin);
c = abs(mean(obj.signal)).^2;
s = mean(abs(obj.signal).^2);
cspr_ = 10*log10(c / s);
cspr =carrier_power_dbm-signal_power_dbm;
end
function papr = papr(obj)
average_power = obj.power;
peak_power = pow2db(max(abs(obj.signal.^2)))+30;
papr = peak_power - average_power;
end
end
end

577
Classes/00_signals/Signal.m
View File

@@ -5,14 +5,21 @@ classdef Signal
properties
signal
logbook
fs
end
methods
function obj = Signal(signal)
function obj = Signal(signal,options)
%SIGNAL Construct an instance of this class
% Detailed explanation goes here
arguments
signal
options.fs = [];
end
obj.signal = signal;
obj.signal = obj.signal;
obj.fs = options.fs;
SignalType = [];
TimeStamp = [];
@@ -114,6 +121,64 @@ classdef Signal
end
function plot(obj, options)
% signal to plot: obj.signal
% fsamp : obj.fs (e.g. 92e9 => 92 GHz)
% length: length(obj.signal)
arguments
obj
options.fignum
options.displayname = [];
options.timeframe = 0;
end
figure(options.fignum); % If figure does not exist, create new figure
% 2) Plot into the figure handle found or created in one
t = (0:length(obj.signal)-1) / obj.fs; % time vector
if options.timeframe ~= 0
%only show a certain timeframe of signal
t = t(t<options.timeframe);
end
% 2 a) Actual plot (hold on, displayname??)
dn = options.displayname;
if isa(obj,'Opticalsignal')
sig = abs(obj.signal).^2;
else
sig = obj.signal;
end
hold on;
plot(t, sig(1:length(t)), 'DisplayName', dn, 'LineWidth', 0.1, 'Marker', '.', 'LineStyle','none', 'MarkerSize', 0.1);
% 2 c)
% - xlabel if not already here: time in readable format (1 ms and not 1e-3 s)
% - ylabel amplitude
if isempty(get(gca, 'XLabel').String)
xlabel('Time (mu s)');
end
if isempty(get(gca, 'YLabel').String)
ylabel('Amplitude');
end
% Convert time axis to milliseconds for readability
xticks = get(gca, 'XTick');
set(gca, 'XTick', xticks, 'XTickLabel', xticks * 1e6);
% Add legend if not already present
if isempty(get(gca, 'Legend'))
legend;
end
hold off;
end
%% Add signals from one signal to another, the first object will sustain
function Sum = plus(X,y)
@@ -126,9 +191,9 @@ classdef Signal
end
end
function Product = times(X,y)
if isa(y,'Signal')
Product = X;
Product.signal = X.signal .* y.signal;
@@ -138,7 +203,6 @@ classdef Signal
end
end
%% Display length
function return_length = length(obj)
@@ -162,7 +226,7 @@ classdef Signal
SignalPower = [obj.power];
Nase = [0];
cell = {SignalType , TimeStamp , Length , SignalPower , Nase, Description};
cell = {SignalType , TimeStamp , Length , SignalPower(1) , Nase, Description};
obj.logbook = [obj.logbook;cell];
@@ -175,9 +239,11 @@ classdef Signal
obj Signal
options.fs_in double
options.fs_out double
options.n double = 10;
options.beta double = 5;
end
obj.signal = resample(obj.signal,options.fs_out,options.fs_in);
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' ];
@@ -188,108 +254,110 @@ classdef Signal
end
%%
function spectrum(obj,fsamp,options)
function spectrum(obj,options)
arguments
obj
fsamp
options.figurename = [];
options.displayname = [];
options.fignum
options.displayname = "";
end
%Get figure if there is already a spectrum plot -> I want to add the new
%spectum "onto" the existing plot to have a better comparison
if isempty(options.figurename)
fig = findall(groot, 'Type', 'figure', 'Name', 'power density');
if isvalid(fig)
fig = get(fig);
ax = gca;
hold on
else
figure('name','power density');
ax = gca;
end
else
fig = findall(groot, 'Type', 'figure', 'Name', options.figurename);
if isvalid(fig)
ax = fig.CurrentAxes;
hold on
else
figure('name',options.figurename);
ax = gca;
hold on
end
end
% spectrum_plot(obj.signal,options.fsamp,options.figurename,options.displayname);
N = 2^(nextpow2(length(obj.signal))-6);
[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"
%compute FFT of input
Fsignal = fft(obj.signal);
%POWER spectral density (todo: toggle?)
psd = Fsignal.*conj(Fsignal);
%Use only magnitude of FFT (which was complex)
psd = abs(psd);
%Shift the spectrum to yield
psd = fftshift(psd);
%divide by N
psd = psd/length(Fsignal);
%smoothing
psd = smooth(psd,1000);
psd_plot = 20*log10(psd);
psd_plot(psd_plot<-120) = -120;
testParseval = 1;
if testParseval == 1
E_FreqDomain = sum(psd);
%test parseval
E_TimeDomain = sum(abs(Fsignal.^2));
if isequal(round(E_FreqDomain,1),round(E_TimeDomain,1))
%disp('Parseval is right!');
else
disp('Parseval theorem is not right...');
end
end
if fsamp <= 1e+100
%Frequency Axis
freq_vec = linspace(-fsamp/2,fsamp/2,length(psd));
freq_vec = reshape(freq_vec,size(psd_plot));
if ~isempty(options.displayname)
plot(freq_vec*1e-9,psd_plot,'Linewidth',0.5,'DisplayName',options.displayname,'Parent',ax);
else
plot(freq_vec*1e-9,psd_plot,'Linewidth',0.5,'Parent',ax);
end
xlabel('Frequency [GHz]')
else
%Wavelength Axis
freq_vec = physconst('LightSpeed')*linspace(-fsamp/2,fsamp/2,length(psd))./((physconst('LightSpeed')/1550e-9)^2);
if ~isempty(options.displayname)
plot(freq_vec*1e9,psd_plot,'Linewidth',0.5,'DisplayName',options.displayname,'Parent',ax);
else
plot(freq_vec*1e9,psd_plot,'Linewidth',0.5,'Parent',ax);
end
xlabel('Wavelength [nm]')
end
ylabel('Magnitude [dB]')
figure(options.fignum); % If figure does not exist, create new figure
hold on
plot(w.*1e-9,p_dbm,'DisplayName',options.displayname,'LineWidth',1);
xlabel("Frequency in GHz");
%ylabel("Power/frequency (dB/Hz)");
ylabel("Power (dBm)");
xlim([-obj.fs/2 obj.fs/2].*1e-9)
edgetick = 2^(nextpow2(obj.fs*1e-9));
xticks([-edgetick:16:edgetick]);
xlim([-244, 244])
ylim([-120,-0]);
yticks([-200:10:10]);
legend
grid minor;
end
%% Power of signal
function pow = power(obj,options)
arguments
obj
options.unit power_notation = power_notation.dBm
end
pow = mean(abs(obj.signal).^2);
switch options.unit
case power_notation.dBm
if isa(obj,'Electricalsignal')
pow = pow / 50;
end
pow = 10*log10(pow)+30; %dbm
case power_notation.mW
pow = pow .* 1e3; %mW
case power_notation.W
%pow = pow % Watt
end
end
%% Peak Power of Signal
function pow_pk = power_peak(obj,options)
arguments
obj
options.unit power_notation = power_notation.dBm
end
pow_pk = max(abs(obj.signal).^2); % dBm
switch options.unit
case power_notation.dBm
pow_pk = pow2db(pow_pk)+30; %dbm
case power_notation.mW
pow_pk = pow_pk .* 1e3; %mW
case power_notation.W
%pow = pow % Watt
end
end
%% PAPR of signal
function papr = papr_lin(obj)
%PAPR The peak-to-average power ratio (PAPR) is the peak amplitude squared (giving the peak power)
% divided by the RMS value squared (giving the average power).[1] It is the square of the crest factor.
% papr = max(abs(timesignal))^2 / rms(timesignal)^2; ODER papr = peak2rms(sig)^2;
papr = obj.power_peak("unit",power_notation.W) / obj.power("unit",power_notation.W); %linear
end
%% PAPR of signal
function papr_db = papr_db(obj)
%PAPR The peak-to-average power ratio (PAPR) is the peak amplitude squared (giving the peak power)
% divided by the RMS value squared (giving the average power).[1] It is the square of the crest factor.
% papr = max(abs(timesignal))^2 / rms(timesignal)^2; ODER papr = peak2rms(sig)^2;
papr = obj.power_peak("unit",power_notation.W) / obj.power("unit",power_notation.W);
papr_db = 10*log10(papr);
average_power = obj.power;
peak_power = obj.power_peak;
papr_db = peak_power - average_power; %db
end
%%
function obj = normalize(obj,options)
@@ -366,33 +434,310 @@ classdef Signal
end
function eye(obj,fsym)
function er = extinctionratio(obj,fsym,M)
histpoints = 1024; %% verticale resolution
histpoints = floor(histpoints/2)*2+1; %% to have the eye digram centered around one point make the vertical resolution uneven
histpoints_horizontal = 512; %% horizontal resolution
hist_data=zeros(histpoints,histpoints_horizontal ); %% initilize eye diagram
disp('h');
fsig = obj.fs;
q = fsig/fsym;
if q > 10 && isinteger(q)
sig = (obj.signal);
if isa(obj,'Opticalsignal')
sig = abs(obj.signal).^2;
elseif isa(obj,'Electricalsignal')
sig = obj.signal;
else
sig = (obj.resample("fs_in",fsig,"fs_out",fsym*30).signal);
q = 10;
sig = obj.signal;
end
figure()
x = (sig); %% make input signal rea)l
x = resample(x,fsym*histpoints_horizontal/2,obj.fs); %% up sample to original fsym rate
if mod(length(x),2)==1 %% if the signal lenght is not divisible by 2 (symbols displayed in the eye diagram are 2) remove last symbol
x = x(1:end-1);
end
eye_mat = reshape(x(1:end-mod(length(x),histpoints_horizontal)),histpoints_horizontal,floor(length(x)/histpoints_horizontal)); %% reshape signal into 256 rows each row has the histogram(eye data of all symbols)
maxA = max(sig(100:end-100));
minA = min(sig(100:end-100));
difference= maxA-minA;
data_ind_y=round((eye_mat-minA)/difference*(histpoints-1)) +1;
for n=1:size(data_ind_y,1)
nn=histcounts(data_ind_y(n,:),1:histpoints+1);
hist_data(:,n)=flip(nn.'); %without flip, the eye is upside down :-(
end
plot_data = 20*log10(hist_data);
plot_data(plot_data==-Inf) = 0;
maxall = 0;
for l = 1:size(plot_data,2)
[maxpk_,pos_] = max(plot_data(:,l));
if maxpk_ > maxall
maxall = maxpk_;
posxall = l;
posyall = pos_;
end
end
hist_interest = plot_data(:,posxall);
hist_interest_smoth = smooth(hist_interest,20);
[pk,loc] = findpeaks(hist_interest_smoth,"MinPeakDistance",40,"NPeaks",M,"MinPeakHeight",30);
for i = 1:numel(loc)
ppeak(i) = maxA - (difference/histpoints*loc(i));
end
if isa(obj,'Opticalsignal')
er=10*log10(ppeak(1)/ppeak(end));
elseif isa(obj,'Electricalsignal')
if mean([ppeak(1),ppeak(end)]) < 1e-2
disp("No Extiction Ration for Bipolar Electrical Signal. Calculating Outer OMA instead...")
er=max(ppeak)-min(ppeak);
else
er=10*log10(ppeak(1)/ppeak(end));
end
else
er=10*log10(ppeak(1)/ppeak(end));
end
if 0
findpeaks(hist_interest_smoth,"MinPeakDistance",40,"NPeaks",M,"MinPeakHeight",30);
end
end
function eye(obj,fsym,M)
mode = 1;
histpoints = 1024; %% verticale resolution
histpoints = floor(histpoints/2)*2+1; %% to have the eye digram centered around one point make the vertical resolution uneven
histpoints_horizontal = 512; %% horizontal resolution
hist_data=zeros(histpoints,histpoints_horizontal ); %% initilize eye diagram
if isa(obj,'Opticalsignal')
sig = abs(obj.signal).^2;
elseif isa(obj,'Electricalsignal')
sig = obj.signal;
else
sig = obj.signal;
end
x = (sig); %% make input signal rea)l
x = resample(x,fsym*histpoints_horizontal/2,obj.fs); %% up sample to original fsym rate
if mod(length(x),2)==1 %% if the signal lenght is not divisible by 2 (symbols displayed in the eye diagram are 2) remove last symbol
x = x(1:end-1);
end
eye_mat = reshape(x(1:end-mod(length(x),histpoints_horizontal)),histpoints_horizontal,floor(length(x)/histpoints_horizontal)); %% reshape signal into 256 rows each row has the histogram(eye data of all symbols)
figure(922)
clf
cursor = 200*q;
if mode == 2
% generate "intuitive eye diagram" by drawing lines on top over
% each other; only draw 1000 lines, otherwise the plot is too
% crowded
col = cbrewer2('Set1',2);
for n=1:1000
hold on
plot(eye_mat(:,n),'LineStyle',':','LineWidth',0.1,'Color',col(2,:));
end
ylabel('Amplitude of Signal');
elseif mode == 1
% generate eye diagram using histogram
maxA = max(sig(100:end-100));
minA = min(sig(100:end-100));
difference= maxA-minA;
data_ind_y=round((eye_mat-minA)/difference*(histpoints-1)) +1;
for n=1:size(data_ind_y,1)
nn=histcounts(data_ind_y(n,:),1:histpoints+1);
hist_data(:,n)=flip(nn.'); %without flip, the eye is upside down :-(
end
plot_data = 20*log10(hist_data);
plot_data(plot_data==-Inf) = 0;
imagesc(plot_data);
% beautify
colormap(cbrewer2("Blues",4096));
if isa(obj,'Opticalsignal')
title("Optical Eye")
ylabel("Power in mW");
y_tickstring = string(linspace(maxA.*1e3,minA.*1e3,16));
min_ = min(abs(obj.signal(100:end-100)).^2);
max_ = abs(max(obj.signal(100:end-100)).^2);
elseif isa(obj,'Electricalsignal')
title("Electrical Eye")
ylabel("Voltage in V");
y_tickstring = string(linspace(maxA,minA,16));
min_ = min(obj.signal(100:end-100));
max_ = abs(max(obj.signal(100:end-100)));
else
title("Digital Eye")
ylabel("Digital Signal Amplitude");
y_tickstring = string(linspace(maxA,minA,16));
min_ = min(obj.signal(100:end-100));
max_ = abs(max(obj.signal(100:end-100)));
end
% add information
if 1
pwr_dbm = round(obj.power,3);
pwr_lin = obj.power("unit",power_notation.W);
papr_ = obj.papr_lin;%round(papr(obj.signal.^2),3);
yline( histpoints-(pwr_lin - minA)/difference*histpoints );
yline( histpoints-(min_ - minA)/difference*histpoints );
yline( histpoints-(max_ - minA)/difference*histpoints );
maxall = 0;
for l = 1:size(plot_data,2)
[maxpk_,pos_] = max(plot_data(:,l));
if maxpk_ > maxall
maxall = maxpk_;
posxall = l;
posyall = pos_;
end
end
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]
boxColor = [0.9 0.9 0.9]; % Light grey background color
boxEdgeColor = 'k'; % Black edge color
boxLineStyle = '--'; % Dashed line style
boxFontWeight = 'bold'; % Bold font
% Create first annotation box for Power
annotation('textbox', boxPosition, ...
'String', ['Power: ',num2str(pwr_dbm),' dBm'], ...
'BackgroundColor', boxColor, ...
'EdgeColor', boxEdgeColor, ...
'LineStyle', boxLineStyle, ...
'FontWeight', boxFontWeight, ...
'HorizontalAlignment', 'center');
% Adjust position for the second box (slightly to the right)
boxPosition = [0.37 0.86 0.2 0.05]; % Adjusted position
% Create second annotation box for PAPR
annotation('textbox', boxPosition, ...
'String', ['PAPR(lin):',num2str(papr_),''], ...
'BackgroundColor', boxColor, ...
'EdgeColor', boxEdgeColor, ...
'LineStyle', boxLineStyle, ...
'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
% 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');
yticks(linspace(0,histpoints,16));
yticklabels(y_tickstring);
end
for s = 1:700
plot(sig(cursor-q:cursor+q),'Color','black','LineWidth',0.1,'LineStyle','-');
hold on
cursor=cursor+q;
s=s+1;
end
ylim([-3 3]);
% disp('h');
%
% fsig = obj.fs;
% q = fsig/fsym;
%
% if q > 10 && isinteger(q)
% sig = (obj.signal);
% else
% sig = (obj.resample("fs_in",fsig,"fs_out",fsym*30).signal);
% q = 10;
% end
%
% figure()
% clf
% cursor = 200*q;
%
% for s = 1:700
% plot(sig(cursor-q:cursor+q),'Color','black','LineWidth',0.1,'LineStyle','-');
% hold on
% cursor=cursor+q;
% s=s+1;
% end
%
% ylim([-3 3]);
end

96
Classes/01_transmit/AWG.m
View File

@@ -46,7 +46,7 @@ classdef AWG < handle
options.lpf_type = 0;
options.f_cutoff = 32e9;
options.H_lpf Filter
end
fn = fieldnames(options);
@@ -60,26 +60,33 @@ classdef AWG < handle
function signalclass_out = process(obj,signalclass_in)
% 0) START AWG SIMULATION
len_in = length(signalclass_in.signal);
% disp(['Power Raw Input: ', num2str(mean(abs(real(signalclass_in.signal).^2)))]);
% 1) RESAMPLE IF NESSECARY
if signalclass_in.fs ~= obj.fdac
k = obj.kover*obj.fdac / signalclass_in.fs;
signalclass_in = signalclass_in.resample("fs_in",signalclass_in.fs,"fs_out",obj.fdac);
% min_ = min(signalclass_in.signal);
% max_ = max(signalclass_in.signal);
signalclass_in = signalclass_in.resample("fs_in",signalclass_in.fs,"fs_out",obj.fdac,"n",10,"beta",5);
% signalclass_in.signal = clip(signalclass_in.signal,-1.3414,1.3414);
% signalclass_in.signal = softclip(signalclass_in.signal,7,1.3414);
% signalclass_in.plot("displayname",'resampled and clipped','fignum',1);
end
% 1-3. actual processing of the signal (normalize->quantize->sample hold)
% disp(['Power after resamp: ', num2str(mean(abs(real(signalclass_in.signal).^2)))]);
% 4) PROCESSING of the signal (normalize->quantize->sample hold)
signalclass_in.signal = obj.process_(signalclass_in.signal);
% cast the inform. signal to electrical signal
if ~isa(signalclass_in,'Electricalsignal')
signalclass_in = Electricalsignal(signalclass_in,"fs",obj.fdac*obj.kover,"logbook",signalclass_in.logbook);
end
% normalize to 0dBm before applying the lowpass
%signalclass_in = signalclass_in.normalize("mode","milliwatt");
signalclass_in = signalclass_in.setPower(13,"dBm");
% 4. Apply LPF on the signal
if obj.lpf_active
if isa(obj.H_lpf,'Filter')
@@ -92,9 +99,9 @@ classdef AWG < handle
signalclass_in = lpf.process(signalclass_in);
end
end
disp(['AWG output power: ',num2str(signalclass_in.power),' dBm']);
% disp(['AWG output power: ',num2str(signalclass_in.power),' dBm']);
% append to logbook
current_class = class(obj);
lbdesc = ['AWG ', current_class , '// k_over:',num2str(obj.kover),'. f_dac:',num2str(obj.fdac*1e-9),'GHz. Resolution:',num2str(obj.bit_resolution),' bits.'];
@@ -125,9 +132,57 @@ classdef AWG < handle
obj.signal_length = length(data_in);
%%%%%%%%% PRECOMP SINC ROLLOFF %%%%%%%%%
if 1
% X: design FIR filter for sinc precomp
ntaps = 13;
npts = 32;
% least-squares FIR design
fmax = obj.fdac*0.5;
ff = linspace(0,fmax,npts);
hsinc = sin(pi*ff/obj.fdac)./(pi*ff/obj.fdac + eps); % transfer function of sample and hold DAC
hsinc(1) = 1;
h_goal= 1./hsinc; % goal function
f = 2.*ff./obj.fdac; %vector between 0 and 1, where 1 is nyquist is fsamp/2
b = firls(ntaps-1,f,h_goal);
data_in = conv(data_in,b,"same");
end
if 0
%compare different fir construction methods in matlab
b_fir = fir2(ntaps-1,f,h_goal);
b_pm = firpm(ntaps-1,f,h_goal);
b = firls(ntaps-1,f,h_goal);
figure(111)
freqz(b,1,[],obj.fdac);
hold on
freqz(b_fir,1,[],obj.fdac);
hold on
freqz(b_pm,1,[],obj.fdac);
plot(f.*obj.fdac./2.*1e-9,20*log10(h_goal),'Marker','o');
legend('firls','fir2','firpm','target')
end
if 0
% design filter as inverse of the goal transfer function
freq_vec = linspace(-obj.fdac/2,obj.fdac/2,length(data_in));
h = sin(pi*freq_vec/obj.fdac)./(pi*freq_vec/obj.fdac + eps);
max_amp_db = 45;
max_amp_lin = 10^(-max_amp_db/20);
p = find(h<max_amp_lin);
% h(p)=10^(-max_amp_db/20);
h_sinc = 1./h;
convol = fftshift(h_sinc).'.*fft(data_in);
data_in = ifft(convol);
data_in = real(data_in);
end
%%%%%%%%% Normalize to DAC range %%%%%%%%%
% X:
if obj.normalize2dac
% 0a Normalize the signal to full scale DAC range
data_in = data_in - min(data_in); %set "foot" to zero
data_in = data_in / (max(data_in)-min(data_in)); %scale between 0 and 1
data_in = data_in * (obj.dac_max-obj.dac_min); %scale to desired total range of DAC
@@ -136,7 +191,10 @@ classdef AWG < handle
data_in = data_in - mean(data_in);
% 1. Quantize the signal - Full Scale is between obj.dac_min
% disp(['Power after scaling to DAC: ', num2str(mean(abs(real(data_in).^2)))]);
%%%%%%%%% Quantize %%%%%%%%%
% X. Quantize the signal - Full Scale is between obj.dac_min
% and dac_max. If signal is smaller in between, you won't use
% the full bit-resolution.
if obj.bit_resolution>0
@@ -144,9 +202,9 @@ classdef AWG < handle
else
elec_out = data_in;
end
% 2. Sample and hold + repeat (data_out: 1xsignal length)
%%%%%%%%% Sample and hold %%%%%%%%%
% X. Sample and hold + repeat (data_out: 1xsignal length)
if obj.upsampling_method == 1
% just use matlab function
elec_out = resample(elec_out,obj.kover,1);
@@ -157,8 +215,8 @@ classdef AWG < handle
else
error('chosen upsampling method not implemented?');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 3. Add skew (not implemented so far)
if obj.skew_active
elec_out = obj.skew(elec_out);

22
Classes/01_transmit/M8196A.m
View File

@@ -5,26 +5,22 @@ classdef M8196A < AWG
methods
function obj = M8196A()
function obj = M8196A(options)
%This is a ready to use AWG simulation that resembles the
%properties of the Keysighe M8196A % M8196A (92GBd) https://www.keysight.com/us/en/product/M8196A/92-gsa-s-arbitrary-waveform-generators.html
arguments
%options.bla = 1;
options.kover = 8;
end
obj = obj@AWG();
obj.dac_max = 0.5;
obj.dac_min = -.5;
obj.bit_resolution = 5.5;
dac_max = 0.4;
dac_min = -0.4;
obj.fdac = 92e9;
obj.lpf_active = 1;
obj.f_cutoff = 32e9;
obj.lpf_type = filtertypes.gaussian;
obj.H_lpf = Filter('filtdegree',5,"f_cutoff",obj.f_cutoff,"fsamp",obj.kover*obj.fdac,"filterType",obj.lpf_type);
fdac = 92e9;
Lp_awg = Filter('filtdegree',4,"f_cutoff",32e9,"fs",fdac*options.kover,"filterType",filtertypes.butterworth,"active",true);
obj = obj@AWG("fdac",fdac,"dac_min",dac_min,"dac_max",dac_max,"lpf_active",1,"H_lpf",Lp_awg,"kover",options.kover,...
"bit_resolution",5.5,"normalize2dac",1,"upsampling_method","samplehold");
end

32
Classes/01_transmit/M8199A.m Normal file
View File

@@ -0,0 +1,32 @@
classdef M8199A < AWG
properties
end
methods
function obj = M8199A(options)
%This is a ready to use AWG simulation that resembles the
%properties of the Keysight M8199A
arguments
options.kover = 8;
end
dac_max = 0.3;
dac_min = -0.3;
fdac = 256e9;
Lp_awg = Filter('filtdegree',4,"f_cutoff",65e9,"fs",fdac*options.kover,"filterType",filtertypes.butterworth,"active",true);
obj = obj@AWG("fdac",fdac,"dac_min",dac_min,"dac_max",dac_max,"lpf_active",1,"H_lpf",Lp_awg,"kover",options.kover,...
"bit_resolution",5.5,"normalize2dac",1,"upsampling_method","samplehold");
end
end
end

21
Classes/01_transmit/M8199B.m
View File

@@ -5,25 +5,22 @@ classdef M8199B < AWG
methods
function obj = M8199B()
function obj = M8199B(options)
%This is a ready to use AWG simulation that resembles the
%properties of the Keysighe M8199B https://www.keysight.com/us/en/assets/3120-1465/data-sheets/M8199A-128-256-GSa-s-Arbitrary-Waveform-Generator.pdf
arguments
%options.bla = 1;
options.kover = 8;
end
obj = obj@AWG();
obj.dac_max = 0.5;
obj.dac_min = -.5;
obj.bit_resolution = 5.5;
dac_max = 0.6;
dac_min = -0.6;
obj.fdac = 256e9;
fdac = 256e9;
Lp_awg = Filter('filtdegree',4,"f_cutoff",75e9,"fs",fdac*options.kover,"filterType",filtertypes.butterworth,"active",true);
obj.lpf_active = 1;
obj.f_cutoff = 80e9;
obj.lpf_type = filtertypes.gaussian;
obj.H_lpf = Filter('filtdegree',5,"f_cutoff",obj.f_cutoff,"fsamp",obj.kover*obj.fdac,"filterType",obj.lpf_type);
obj = obj@AWG("fdac",fdac,"dac_min",dac_min,"dac_max",dac_max,"lpf_active",1,"H_lpf",Lp_awg,"kover",options.kover,...
"bit_resolution",5.5,"normalize2dac",1,"upsampling_method","samplehold");
end

117
Classes/01_transmit/PAMsource.m Normal file
View File

@@ -0,0 +1,117 @@
classdef PAMsource
%NAME Summary of this class goes here
% Detailed explanation goes here
properties(Access=public)
order
useprbs
M
fsym
randkey
applypulseform
pulseformer
fs_out
applyclipping
clipfactor
end
methods (Access=public)
function obj = PAMsource(options)
%NAME Construct an instance of this class
% Detailed explanation goes here
arguments
options.order = 16;
options.useprbs = true;
options.M = true
options.fsym = 112e9;
options.randkey = 0;
options.applypulseform = 1;
options.pulseformer Pulseformer;
options.fs_out ;
options.applyclipping = 0;
options.clipfactor = 10;
end
fn = fieldnames(options);
for n = 1:numel(fn)
try
obj.(fn{n}) = options.(fn{n});
end
end
if boolean(obj.applypulseform) && isempty(obj.pulseformer)
warning('No Pulseformer given. Proceeding with RRC and alpha 0.05');
options.pulseformer = Pulseformer("fsym",obj.fsym,"fdac",obj.fs_out,"pulse","rrc","pulselength",16,"rrcalpha",0.05);
end
end
function [digi_sig,symbols,bits] = process(obj)
%%%%% PRBS Generation in correct shape for Modulation Format %%%%%%
O = obj.order; %order of prbs
N = 2^(O-1); %length of prbs
[~,seed] = prbs(O,1); %initialize first seed of prbs
bitpattern=[];
if obj.useprbs
for i = 1:log2(obj.M)
[bitpattern(:,i),seed] = prbs(O,N,seed);
end
else
s = RandStream('twister','Seed',obj.randkey);
for i = 1:log2(obj.M)
bitpattern(:,i) = randi(s,[0 1], N, 1);
end
end
if obj.M == 6
bitpattern = reshape(bitpattern,[],1);
bitpattern = bitpattern(1:end-mod(length(bitpattern),5));
end
bits = Informationsignal(bitpattern);
symbols = PAMmapper(obj.M,0).map(bits);
symbols.fs = obj.fsym;
if obj.applyclipping
sym_min = min(symbols.signal);
sym_max = max(symbols.signal);
end
%%%%% Pulseforming %%%%%%
if obj.applypulseform
digi_sig = obj.pulseformer.process(symbols);
else
digi_sig = symbols;
end
%%%%% Resample to f DAC %%%%%%
digi_sig = digi_sig.resample("fs_in",digi_sig.fs,"fs_out",obj.fs_out,"n",10,"beta",5);
%%%%% Hard clip digital signal to PAM range before DAC %%%%%%
if obj.applyclipping
digi_sig.signal = clip(digi_sig.signal , sym_min * obj.clipfactor , sym_max * obj.clipfactor);
end
% append to logbook
lbdesc = ['Generated PAM Signal'];
digi_sig = digi_sig.logbookentry(lbdesc);
end
end
end

10
Classes/01_transmit/Pulseformer.m
View File

@@ -123,9 +123,9 @@ classdef Pulseformer
%data_out_ = upfirdn(data_in,h,up,dn);
%
% %cut signal, which is longer due to fir filter
st = round(up/dn*racos_len/2); %we need to cut y_out
% en = round(st + (length(data_in)*up/dn) -1);
% data_out = data_out(st:en);
% st = round(up/dn*racos_len/2); %we need to cut y_out
% en = round(st + (length(data_in)*up/dn) -1);
% data_out = data_out(st:en);
%scaling?! see pulsef module line 696
% scale = max(max([abs(real(data_out)) abs(imag(data_out))])); %find max value from real and imag part
@@ -134,8 +134,8 @@ classdef Pulseformer
data_out = data_out';
%Check output integrity
if round(up/dn * length(data_in)) ~= length(data_out)
%warning('Check signal length after pulse shaping');
if abs(round(up/dn * length(data_in)) - length(data_out)) > 4
warning('Check signal length after pulse shaping');
%disp('Check signal length after pulse shaping');
end

87
Classes/01_transmit/Signalgenerator.m Normal file
View File

@@ -0,0 +1,87 @@
classdef Signalgenerator
%NAME Summary of this class goes here
% Detailed explanation goes here
properties(Access=public)
form
length
fs
fsig
end
methods (Access=public)
function obj = Signalgenerator(options)
%NAME Construct an instance of this class
% Detailed explanation goes here
arguments
options.form signalform = signalform.sine
options.length double = 1024
options.fs double = 1000 %Hz sampling
options.fsig double = 50 % Hz fundamental frex e.g. of the sine or sawtooth
end
%
fn = fieldnames(options);
for n = 1:numel(fn)
try
obj.(fn{n}) = options.(fn{n});
end
end
end
function signalclass_out = process(obj)
% actual processing of the signal (steps 1. - 3.)
signal = obj.build_signal();
signalclass_out = Informationsignal(signal,"fs",obj.fs);
% append to logbook
lbdesc = ['Signalgenerator: ',char(obj.form), ''];
signalclass_out = signalclass_out.logbookentry(lbdesc);
end
function signal = build_signal(obj)
%METHOD1 Summary of this method goes here
% Detailed explanation goes here
arguments(Input)
obj
end
arguments(Output)
signal double
end
switch obj.form
case signalform.sine
% Parameters
T = 1/obj.fs; % Sampling period (seconds per sample)
L = obj.length; % Length of signal (number of samples)
t = (0:L-1)*T; % Time vector
% Sine wave parameters
f = obj.fsig; % Frequency of the sine wave (Hz)
A = 1; % Amplitude of the sine wave
% Generate sine wave
signal = A * sin(2*pi*f*t);
case signalform.noise
end
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

12
Classes/02_etc/Filter.m
View File

@@ -3,6 +3,8 @@ classdef Filter < handle
% Detailed explanation goes here
properties
active
H
filterType
f_cutoff
@@ -22,6 +24,7 @@ classdef Filter < handle
%FILTER Construct an instance of this class
% Detailed explanation goes here
arguments
options.active logical = true;
options.filterType filtertypes = filtertypes.bessel_inp ;
options.f_cutoff = 0;
options.fs = 0;
@@ -47,7 +50,10 @@ classdef Filter < handle
end
function signal_out = process(obj,signal_in)
if ~obj.active
signal_out = signal_in;
return
end
if isa(signal_in,'Signal')
% actual processing of the signal
@@ -284,10 +290,12 @@ classdef Filter < handle
title('Magnitude')
hold on
p = plot(frex,filt,'LineWidth',1,'LineStyle','-','DisplayName',[cur_module_name,': ',filter_desc, '; 3/6/10 dB: ', num2str(threedB),'/ ',num2str(sixdB),'/ ',num2str(ninedB)]);
if 0
xline([-threedB,threedB],'LineStyle','-','LineWidth',1,'HandleVisibility','off','Color',p.Color);
xline([-sixdB,sixdB],'LineStyle','-','LineWidth',1,'HandleVisibility','off','Color',p.Color);
xline([-ninedB,ninedB],'LineStyle','-','LineWidth',1,'HandleVisibility','off','Color',p.Color);
end
xline(fcut_,'LineStyle',':','LineWidth',1,'HandleVisibility','off','Color',p.Color);
yline([-3, -6, -9],'LineStyle',':','LineWidth',1,'HandleVisibility','off');

42
Classes/02_optical/Fiber.m
View File

@@ -49,7 +49,7 @@ classdef Fiber
function signalclass_out = process(obj,signalclass_in)
% actual processing of the signal (steps 1. - 3.)
signalclass_in.signal = obj.process_(signalclass_in.signal);
signalclass_in = obj.process_(signalclass_in);
% append to logbook
lbdesc = 'Fiber ';
@@ -60,27 +60,45 @@ classdef Fiber
end
function opt_out = process_(obj,opt_in)
function opt_sig = process_(obj,opt_sig)
%METHOD1 Summary of this method goes here
% Detailed explanation goes here
obj.b2 = -obj.D*obj.lambda0^2/(2*pi*Constant.LightSpeed);
obj.b3 = ((obj.lambda0.^2/(2*pi*Constant.LightSpeed)).^2*obj.Dslope);
obj.alpha_lin = obj.alpha/10*log(10)/1000;
signal = opt_sig.signal;
lambda_signal = opt_sig.lambda;
N = length(opt_in);
obj.D = obj.D + (lambda_signal-obj.lambda0)*obj.Dslope;
obj.b2 = -obj.D*lambda_signal^2/(2*pi*Constant.LightSpeed);
obj.b3 = ((lambda_signal.^2/(2*pi*Constant.LightSpeed)).^2*obj.Dslope);
obj.alpha_lin = obj.alpha/10*log(10)/1000;
N = length(signal);
faxis = linspace(-obj.fsimu/2,obj.fsimu/2,N+1);
faxis = ifftshift(faxis(:,1:end-1));
faxis = faxis';
obj.linstep = -obj.alpha_lin/2 - 2*1j*pi^2*obj.b2*faxis.^2 - 4/3*1j*pi^3*obj.b3*faxis.^3;
if obj.gamma ~= 0
opt_out = obj.NLSE(opt_in);
else
opt_out = ifft(fft(opt_in).*exp(obj.linstep*obj.fiber_length)); % only one linear step
if 0
H = exp((obj.linstep)*obj.fiber_length);
figure(222)
hold on
plot(faxis.*1e-9,abs(real(Y)),'LineStyle','-','DisplayName','Abs(real) part of complex TF');
xlabel("Frequency in GHz")
ylabel("$ R|(H(\omega, L))|$")
end
%attenuate nase
if obj.gamma ~= 0
opt_out = obj.NLSE(signal);
else
opt_out = ifft( fft(signal) .* exp(obj.linstep*obj.fiber_length) ); % only one linear step
end
%TODO: attenuate nase ...
opt_sig.signal = opt_out;
end

6
Classes/03_receive/Photodiode.m
View File

@@ -41,6 +41,8 @@ classdef Photodiode
% cast the inform. signal to electrical signal
[signalclass_in, nase, lambda] = Electricalsignal(signalclass_in,"fs",obj.fsimu,"logbook",signalclass_in.logbook);
% append to logbook
lbdesc = ['Photo Diode '];
signalclass_in = signalclass_in.logbookentry(lbdesc);
@@ -71,8 +73,8 @@ classdef Photodiode
% Thermal Noise
therm_current_psd = (2 * k * T / R ) ; %squared
Bw = obj.fsimu; %is this correct? shouldnt it be the bandwidth of the actual component? e.g. 70GHz?
therm_noise_pow = therm_current_psd * 2 * Bw; %squared
Bw = obj.fsimu;
therm_noise_pow = therm_current_psd * Bw; %squared
therm_noise = sqrt(therm_noise_pow) .* randn(obj.randomstream,size(yout,1),1);

23
Classes/04_DSP/A1_scheme.m Normal file
View File

@@ -0,0 +1,23 @@
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

192
Classes/04_DSP/EQ_silas.m
View File

@@ -170,23 +170,23 @@ classdef EQ_silas < handle
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;
@@ -194,18 +194,17 @@ classdef EQ_silas < handle
%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
@@ -227,42 +226,29 @@ classdef EQ_silas < handle
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
mu_ffe = [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)];
mu_dfe = ones(1,sum(obj.Cb))*obj.mu_dfe_dd;
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 = zeros(obj.x_length,1);
lvl_err = NaN(obj.x_length,numel(obj.d_constellation));
lvl_err_mov = NaN(100,numel(obj.d_constellation));
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;
if obj.eq_avg_blocklength > 0
averaging_window = zeros(obj.eq_avg_blocklength,1);
end
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).';
%
if obj.eq_avg_blocklength > 0 %% Das läuft gut mit 400er Fenster!!
averaging_window = circshift(averaging_window,obj.sps);
averaging_window(1:obj.sps,1) = x(1:obj.sps);
avg_(k) = mean(averaging_window);
x = x-avg_(k);
end
%bring this signal to "special" VNLE format
x_vnle = obj.calcVNLENonlinVecs(x,obj.Ie2,obj.Ie3,obj.Ne,obj.x_norm);
@@ -270,61 +256,42 @@ classdef EQ_silas < handle
x_d = [x_vnle;-d_vnle];
%Apply filter
if obj.mu_dc_dd > 0
y(m) = obj.e_dc(end) + x_d.'* coeff;
else
y(m) = x_d.'* coeff;
end
y(m) = x_d.'* coeff;
%Decision 1
[~,symbol_idx] = min(abs(y(m) - obj.d_constellation)); % decision for closest constellation point
d_hat(k) = obj.d_constellation(symbol_idx);
%Error between FFE & DFE filtered signal and Decision
obj.error(k) = y(m) - d_hat(k);
%
lvl_err(k,symbol_idx) = obj.error(k);
d_hat(m,symbol_idx) = obj.d_constellation(symbol_idx);
lvl_err_mov(:,symbol_idx) = circshift(lvl_err_mov(:,symbol_idx),1);
lvl_err_mov(1,symbol_idx) = obj.error(k);
y_1(m,symbol_idx) = y(m); % after 1st iteration
%Decision 2
y(m) = y(m)-mean(lvl_err_mov(:,symbol_idx),'omitnan');
[~,symbol_idx] = min(abs(y(m) - obj.d_constellation)); % decision for closest constellation point
d_hat(k) = obj.d_constellation(symbol_idx);
%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(k) * conj(x_d)));
%Update DC error
dc_block(dc_cnt) = obj.error(k) ;
if dc_cnt == obj.eq_parallelization_blocklength
if obj.eq_updatelatency > 1
obj.e_dc = circshift(obj.e_dc,1);
% m_reg(end+1) = ((1:obj.eq_parallelization_blocklength)' \ (cumsum(dc_block)));
%
% obj.e_dc(1) = obj.e_dc(2) - sign(m_reg(end)) .* (sum(dc_block).* m_reg(end) .* obj.mu_dc_dd);
obj.e_dc(1) = obj.e_dc(2) - sum(dc_block) .* obj.mu_dc_dd;
else
%m_reg(end+1) = ((1:obj.eq_parallelization_blocklength)' \ (cumsum(dc_block)));
% obj.e_dc = obj.e_dc - sign(m_reg(end)) .* (sum(dc_block).* m_reg(end) .* obj.mu_dc_dd);
obj.e_dc = obj.e_dc - sum(dc_block) .* obj.mu_dc_dd;
% obj.e_dc = obj.e_dc - obj.mu_dc_dd * obj.error(k); %newapril
end
dc_cnt = 0;
end
coeff = coeff - (mu_mat * (obj.error(m) * conj(x_d)));
% Append new decision to decision feedback
if obj.Nb(1) > 0
@@ -332,43 +299,52 @@ classdef EQ_silas < handle
%shift up one index
d_feedback(2:end) = d_feedback(1:end-1);
%replace 1st index with current estimation
d_feedback(1) = d_hat(k);
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
%%
%
% b = movmean(lvl_err,[500 500],1,"omitnan");
%
% figure(11)
% for i = 1:4
% hold on
% stem(lvl_err(:,i))
% end
%%
obj.y_out = (circshift( y.' ,-(obj.delay))).';
obj.d_out = d_hat(1:2:end);
% err = obj.error(1:2:end);
% res = NaN(8,length(err));
% for lvl = 1:8
% a = find(obj.d_out==obj.d_constellation(lvl));
% res(lvl,a) = err(a);
% end
% mean(res,2,"omitnan");
%
% figure(12)
% scatter(1:length(obj.y_out),obj.y_out,1,'.')
% hold on
% scatter(1:length(obj.d_out),obj.d_out,1,'.')
% 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

456
Classes/04_DSP/EQ_silas_ofc.m Normal file
View File

@@ -0,0 +1,456 @@
classdef EQ_silas_ofc < 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
% 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_ofc(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] = 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';
% append to logbook
lbdesc = ['EQ von Silas ist gelaufen '];
signalclass_in = signalclass_in.logbookentry(lbdesc);
% 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);
lvl_err_true = NaN(length(obj.d),numel(obj.d_constellation));
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
mu_ffe = [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)];
mu_dfe = ones(1,sum(obj.Cb))*obj.mu_dfe_dd;
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 = zeros(obj.x_length,1);
m_reg = 0;
if obj.eq_avg_blocklength > 0
averaging_window = zeros(obj.eq_avg_blocklength,1);
end
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).';
%
if obj.eq_avg_blocklength > 0 %% Das läuft gut mit 400er Fenster!!
averaging_window = circshift(averaging_window,obj.sps);
averaging_window(1:obj.sps,1) = x(1:obj.sps);
avg_(k) = mean(averaging_window);
x = x-avg_(k);
end
%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) = (m_reg(end)*dc_cnt + obj.e_dc(end)) + x_d.'* coeff;
if obj.mu_dc_dd > 0
y(m) = obj.e_dc(end) + x_d.'* coeff;
else
y(m) = x_d.'* coeff;
end
% if obj.eq_avg_blocklength > 0 %% Das läuft nicht gut!!
% averaging_window = circshift(averaging_window,obj.sps);
% averaging_window(1:obj.sps,1) = y(m);
% avg_(m) = mean(averaging_window);
% y(m) = y(m)-avg_(m);
% end
%Decision
[~,symbol_idx] = min(abs(y(m) - obj.d_constellation)); % decision for closest constellation point
d_hat(k) = obj.d_constellation(symbol_idx);
%Error between FFE & DFE filtered signal and Decision
obj.error(k) = y(m) - d_hat(k);
lvl_err_true(m,symbol_idx) = y(m) - obj.d(m+1);
% if obj.eq_avg_blocklength > 0 %% Das läuft nicht gut!!
% averaging_window = circshift(averaging_window,obj.sps);
% averaging_window(1:obj.sps,1) = y(m);
% avg_(m) = mean(averaging_window);
% y(m) = y(m)-avg_(m);
% end
%Update FFE and DFE coefficients
coeff = coeff - mu_mat*obj.error(k) * conj(x_d);
%Update DC error
dc_block(dc_cnt) = obj.error(k) ;
if dc_cnt == obj.eq_parallelization_blocklength
if obj.eq_updatelatency > 1
obj.e_dc = circshift(obj.e_dc,1);
% m_reg(end+1) = ((1:obj.eq_parallelization_blocklength)' \ (cumsum(dc_block)));
%
% obj.e_dc(1) = obj.e_dc(2) - sign(m_reg(end)) .* (sum(dc_block).* m_reg(end) .* obj.mu_dc_dd);
obj.e_dc(1) = obj.e_dc(2) - sum(dc_block) .* obj.mu_dc_dd;
else
%m_reg(end+1) = ((1:obj.eq_parallelization_blocklength)' \ (cumsum(dc_block)));
% obj.e_dc = obj.e_dc - sign(m_reg(end)) .* (sum(dc_block).* m_reg(end) .* obj.mu_dc_dd);
obj.e_dc = obj.e_dc - sum(dc_block) .* obj.mu_dc_dd;
end
dc_cnt = 0;
end
% 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(k);
%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))).';
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)
x1 = x_in_block;
x2 = [];
x3 = [];
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
x2 = input_vec_se(I_2(:,1)).*input_vec_se(I_2(:,2));
end
if N_(3) > 0
delta_3 = round((N_(1)-N_(3))/2);
input_vec_th = x_in_block(delta_3:end)/norm_(3);
x3 = input_vec_th(I_3(:,1)).*input_vec_th(I_3(:,2)).*input_vec_th(I_3(:,3));
end
x_in_vnle_format = [x1;x2;x3];
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

4
Classes/04_DSP/EQ_silas_sliding_window_dc_removal.m
View File

@@ -54,8 +54,6 @@ classdef EQ_silas_sliding_window_dc_removal < handle
eq_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)
end
methods
@@ -110,7 +108,7 @@ classdef EQ_silas_sliding_window_dc_removal < handle
% actual processing of the signal (steps 1. - 3.)
% 1 normalize RMS
%signalclass_in = signalclass_in.normalize("mode","rms");
signalclass_in = signalclass_in.normalize("mode","rms");
% Process the EQ optimization
obj.process_(signalclass_in.signal', reference_signalclass_in.signal');

170
Classes/04_DSP/FFE.m Normal file
View File

@@ -0,0 +1,170 @@
classdef FFE < handle
% Implementation of plain and simple FFE.
% 1) Training mode (stable performance when you use NLMS)
% 2) Decision directed mode
% Eq = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",25,"sps",2,"decide",0);
properties
sps % usually 2
order
e
error
len_tr
mu_tr
epochs_tr
mu_dd
epochs_dd
constellation
decide
end
methods
function obj = FFE(options)
arguments(Input)
options.sps = 2;
options.order = 15;
options.len_tr = 4096;
options.mu_tr = 0;
options.epochs_tr = 5;
options.mu_dd = 1e-5;
options.epochs_dd = 5;
options.decide = false;
end
fn = fieldnames(options);
for n = 1:numel(fn)
obj.(fn{n}) = options.(fn{n});
end
obj.e = zeros(obj.order,1);
obj.error = 0;
end
function [X] = process(obj, X, D)
% actual processing of the signal (steps 1. - 3.)
% 1 normalize RMS
X = X.normalize("mode","rms");
obj.constellation = unique(D.signal);
% Training Mode
training = 1;
showviz = 0;
obj.equalize(X.signal, D.signal,obj.mu_tr,obj.epochs_tr,obj.len_tr,training,showviz);
% Decision Directed Mode
N = X.length;
training = 0;
showviz = 0;
[signal,decision]=obj.equalize(X.signal, D.signal,obj.mu_dd,obj.epochs_dd,N,training,showviz);
% Output Signal
if obj.decide
X.signal = decision;
else
X.signal = signal;
end
X.fs = D.fs; %change sampling frequency of outgoing signal from fdac e.g. 2 sps to symbol spaced = fsym
lbdesc = [num2str(obj.order),' tap FFE'];
X = X.logbookentry(lbdesc); % append to logbook
end
function [y,d_hat] = equalize(obj,x,d,mio,epochs,N,training,showviz)
arguments
obj
x
d
mio
epochs
N
training
showviz
end
x = [zeros(floor(obj.order/2),1); x; zeros(obj.order,1)];
if showviz
f = figure(111);
subplot(2,2,1:2);
hold on
a = scatter(1:numel(x),x,1,'.');
a2 = scatter(1,1,1,'.');
a3 = scatter(1,1,2,'.');
a4 = xline(1);
ylim([-3 3])
xlim([0 length(x)]);
% subplot(2,2,3)
% dplot = x(1:1+500);
% b = scatter(1:numel(dplot),dplot,5,'x');
% xline(1)
% xline(obj.order)
% ylim([-3 3])
% xlim([0 500]);
subplot(2,2,3:4)
c = stem(obj.e);
ylim([-1 1])
drawnow
end
for epoch = 1 : epochs
symbol = 0;
for sample = 1 : obj.sps : N
symbol = symbol+1;
U = x(obj.order+sample-1:-1:sample);
y(symbol,1) = obj.e.' * U; % Calculating output of LMS __ * |
if training
d_hat(symbol,1) = d(symbol);
else
[~,symbol_idx] = min(abs(y(symbol) - obj.constellation)); % decision for closest constellation point
d_hat(symbol,1) = obj.constellation(symbol_idx);
end
err(symbol) = y(symbol) - d_hat(symbol); % Instantaneous error
if mio ~= 0
obj.e = obj.e - (mio * err(symbol) * U) ; % Weight update rule of LMS
else
normalizationfactor = (U.' * U);
obj.e = obj.e - err(symbol) * U / normalizationfactor; % Weight update rule of NLMS
end
if mod(sample,100) == 1 && showviz
a2.XData = 1:2*numel(y);
a2.YData = repelem(y, 2);
a3.XData = 1:2*numel(d_hat);
a3.YData = repelem(d_hat, 2);
a4.Value = sample;
% b.YData = x(symbol:symbol+500);
c.YData = obj.e;
drawnow;
end
obj.error(epoch,symbol) = err(symbol) * err(symbol)'; % Instantaneous square error
end
end
end
end
end

194
Classes/04_DSP/FFE_DFE.m Normal file
View File

@@ -0,0 +1,194 @@
classdef FFE_DFE < handle
% Implementation of plain and simple FFE.
% 1) Training mode (stable performance when you use NLMS)
% 2) Decision directed mode
% Eq = FFE_DFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"ffe_mu_dd",1e-4,"dfe_mu_dd",5e-4,"ffe_mu_tr",0,"dfe_mu_tr",0,"ffe_order",21,"dfe_order",0,"sps",2,"decide",1);
properties
sps % usually 2
ffe_order
dfe_order
e
b
error
len_tr
ffe_mu_tr
dfe_mu_tr
epochs_tr
ffe_mu_dd
dfe_mu_dd
epochs_dd
constellation
decide
end
methods
function obj = FFE_DFE(options)
arguments(Input)
options.sps = 2;
options.ffe_order = 15;
options.dfe_order = 2;
options.len_tr = 4096;
options.ffe_mu_tr = 0;
options.dfe_mu_tr = 0;
options.epochs_tr = 5;
options.ffe_mu_dd = 1e-5;
options.dfe_mu_dd = 1e-5;
options.epochs_dd = 5;
options.decide = false;
end
fn = fieldnames(options);
for n = 1:numel(fn)
obj.(fn{n}) = options.(fn{n});
end
obj.e = zeros(obj.ffe_order,1);
obj.b = zeros(obj.dfe_order,1);
obj.error = 0;
end
function [X] = process(obj, X, D)
% actual processing of the signal (steps 1. - 3.)
% 1 normalize RMS
X = X.normalize("mode","rms");
obj.constellation = unique(D.signal);
% Training Mode
training = 1;
showviz = 0;
obj.equalize(X.signal, D.signal,obj.ffe_mu_tr,obj.dfe_mu_tr,obj.epochs_tr,obj.len_tr,training,showviz);
% Decision Directed Mode
N = X.length;
training = 0;
showviz = 0;
[signal,decision]=obj.equalize(X.signal, D.signal,obj.ffe_mu_dd,obj.dfe_mu_dd,obj.epochs_dd,N,training,showviz);
% Output Signal
if obj.decide
X.signal = decision;
else
X.signal = signal;
end
X.fs = D.fs; %change sampling frequency of outgoing signal from fdac e.g. 2 sps to symbol spaced = fsym
lbdesc = [num2str(obj.ffe_order),' tap FFE'];
X = X.logbookentry(lbdesc); % append to logbook
end
function [y,d_hat] = equalize(obj,x,d,ffe_mu,dfe_mu,epochs,N,training,showviz)
arguments
obj
x
d
ffe_mu
dfe_mu
epochs
N
training
showviz
end
mu = diag([ones(1,obj.ffe_order(1))*ffe_mu(1) ...
ones(1,obj.dfe_order(1))*dfe_mu(1) ]);
x = [zeros(floor(obj.ffe_order/2),1); x; zeros(obj.ffe_order,1)];
%d = [zeros(obj.dfe_order-1,1); d; zeros(obj.dfe_order,1)];
d_ = zeros(obj.dfe_order(1),1);
coeff = [obj.e;obj.b];
if showviz
f = figure(111);
subplot(2,2,1:2);
hold on
a = scatter(1:numel(x),x,1,'.');
a2 = scatter(1,1,1,'.');
a3 = scatter(1,1,2,'.');
a4 = xline(1);
ylim([-3 3])
xlim([0 length(x)]);
subplot(2,2,3:4)
c = stem(obj.e);
ylim([-1 1])
drawnow
end
for epoch = 1 : epochs
symbol = 0;
for sample = 1 : obj.sps : N
symbol = symbol+1;
x_ = x(obj.ffe_order+sample-1:-1:sample);
v = [x_;d_];
y(symbol,1) = coeff.' * v; % Calculating output of LMS __ * |
if training
d_hat(symbol,1) = d(symbol);
else
[~,symbol_idx] = min(abs(y(symbol) - obj.constellation)); % decision for closest constellation point
d_hat(symbol,1) = obj.constellation(symbol_idx);
end
err(symbol) = y(symbol) - d_hat(symbol); % Instantaneous error
if ~all(mu == 0,'all') %not all mu values are zero
coeff = coeff - (mu * err(symbol) * v) ; % Weight update rule of LMS
else
normalizationfactor = (v.' * v);
coeff = coeff - err(symbol) * v / normalizationfactor; % Weight update rule of NLMS
end
% Append new decision to decision feedback
if obj.dfe_order(1) > 0
%shift up one index
d_(2:end) = d_(1:end-1);
%replace 1st index with current estimation
d_(1) = d_hat(symbol);
end
if mod(sample,100) == 1 && showviz
a2.XData = 1:2*numel(y);
a2.YData = repelem(y, 2);
a3.XData = 1:2*numel(d_hat);
a3.YData = repelem(d_hat, 2);
a4.Value = sample;
c.YData = obj.e;
drawnow;
end
obj.error(epoch,symbol) = err(symbol) * err(symbol)'; % Instantaneous square error
end
end
obj.e = coeff(1:obj.ffe_order);
obj.b = coeff(obj.ffe_order+1:end);
end
end
end

293
Classes/04_DSP/VNLE.m Normal file
View File

@@ -0,0 +1,293 @@
classdef VNLE < handle
% Implementation of plain and simple FFE.
% 1) Training mode (stable performance when you use NLMS)
% 2) Decision directed mode
properties
sps % usually 2
order
e
error
len_tr
mu_tr
epochs_tr
mu_dd
epochs_dd
constellation
decide
x_norm
ce
ie2
ie3
end
methods
function obj = VNLE(options)
arguments(Input)
options.sps = 2;
options.order = [15,2,2];
options.len_tr = 4096;
options.mu_tr = 0;
options.epochs_tr = 5;
options.mu_dd = 1e-5;
options.epochs_dd = 5;
options.decide = false;
end
fn = fieldnames(options);
for n = 1:numel(fn)
obj.(fn{n}) = options.(fn{n});
end
obj.error = 0;
end
function [X] = process(obj, X, D)
% actual processing of the signal (steps 1. - 3.)
% 1 normalize RMS
X = X.normalize("mode","rms");
obj.constellation = unique(D.signal);
obj.x_norm = obj.calcPowerNormalization(X.signal);
obj.ce = obj.calcVNLEMemoryLength(obj.order);
[obj.ie2,obj.ie3] = obj.calcIndiceVectors(obj.order);
obj.e = zeros( sum(obj.ce) ,1);
% Training Mode
training = 1;
showviz = 0;
obj.equalize(X.signal, D.signal,obj.mu_tr,obj.epochs_tr,obj.len_tr,training,showviz);
% Decision Directed Mode
N = X.length;
training = 0;
showviz = 0;
[signal,decision]=obj.equalize(X.signal, D.signal,obj.mu_dd,obj.epochs_dd,N,training,showviz);
% Output Signal
if obj.decide
X.signal = decision;
else
X.signal = signal;
end
X.fs = D.fs; %change sampling frequency of outgoing signal from fdac e.g. 2 sps to symbol spaced = fsym
lbdesc = [num2str(obj.order),' tap FFE'];
X = X.logbookentry(lbdesc); % append to logbook
end
function [y,d_hat] = equalize(obj,x,d,mu,epochs,N,training,showviz)
arguments
obj
x
d
mu
epochs
N
training
showviz
end
if all(mu == mu(1))
% mu = mu(1);
mu = diag(ones(1,sum(obj.ce))*mu(1));
else
mu = diag([ones(1,obj.ce(1))*mu(1) ...
ones(1,obj.ce(2))*mu(2) ...
ones(1,obj.ce(3))*mu(3) ]);
end
x = [zeros(floor(obj.order(1)/2),1); x; zeros(obj.order(1),1)];
if showviz
f = figure(111);
subplot(2,2,1:2);
hold on
a = scatter(1:numel(x),x,1,'.');
a2 = scatter(1,1,1,'.');
a3 = scatter(1,1,2,'.');
a4 = xline(1);
ylim([-3 3])
xlim([0 length(x)]);
subplot(2,2,3:4)
c = stem(obj.e);
ylim([-1 1])
drawnow
end
for epoch = 1 : epochs
symbol = 0;
for sample = 1 : obj.sps : N
symbol = symbol+1;
% x_in = x(obj.order(1)+sample+(obj.sps-1):-1:sample+obj.sps);
x_in = x(obj.order(1)+sample-1:-1:sample);
x_in = obj.calcVNLENonlinVecs(x_in,obj.ie2,obj.ie3,obj.order,obj.x_norm);
y(symbol,1) = obj.e.' * x_in; % Calculating output of LMS __ * |
if training
err(symbol) = y(symbol) - d(symbol); % Instantaneous error
else
[~,symbol_idx] = min(abs(y(symbol) - obj.constellation)); % decision for closest constellation point
d_hat(symbol,1) = obj.constellation(symbol_idx);
err(symbol) = y(symbol) - d_hat(symbol); % Instantaneous error
end
if ~all(mu==0,'all') %mu has not only zeros
obj.e = obj.e - (mu * err(symbol) * x_in) ; % Weight update rule of LMS
else
normalizationfactor = (x_in.' * x_in);
obj.e = obj.e - err(symbol) * x_in / normalizationfactor; % Weight update rule of NLMS
end
if mod(sample,100) == 1 && showviz
a2.XData = 1:2*numel(y);
a2.YData = repelem(y, 2);
a3.XData = 1:2*numel(d_hat);
a3.YData = repelem(d_hat, 2);
a4.Value = sample;
% b.YData = x(symbol:symbol+500);
c.YData = obj.e;
drawnow;
end
obj.error(epoch,symbol) = err(symbol) * err(symbol)'; % Instantaneous square error
end
end
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
indvec2nd=[];
indvec3rd=[];
for o = 2:numel(N)
n = N(o);
v = 1:n; % Ursprünglicher Vektor
row = 1;
% Schleifen zur Generierung des Indize Vektors
switch o
case 2
indvec2nd = zeros(n*(n+1)/2, o);
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

BIN
Classes/matlab.mat Normal file
View File

Binary file not shown.