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

small updates from work pc

Browse Source 
try to implement kalman filter for MPI mitigation
This commit is contained in:
sioe
2024-10-15 10:16:35 +02:00
parent 1540f87850
commit d8b4d6fe7c
12 changed files with 799 additions and 121 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
Classes/00_signals/Signal.m
View File

@@ -284,7 +284,7 @@ classdef Signal
N = 2^(nextpow2(length(obj.signal))-8);
[p_lin,w] = pwelch(obj.signal,hanning(N),N/2,N,obj.fs,"centered","power","mean");
normalize = 1;
normalize = 0;
if normalize
p_lin = p_lin./ max(p_lin);
p_dbm = 10*log10(p_lin); %dB to dBm in case of "power"

7
Classes/04_DSP/Equalizer/EQ.m
View File

@@ -104,10 +104,10 @@ classdef EQ
end
function [signalclass_out,error_log] = process(obj,signalclass_in, reference_signalclass_in)
function [signalclass_out,noi] = process(obj,signalclass_in, reference_signalclass_in)
% actual processing of the signal (steps 1. - 3.)
[signalclass_in.signal,error_log] = obj.process_(signalclass_in.signal', reference_signalclass_in.signal');
[signalclass_in.signal] = obj.process_(signalclass_in.signal', reference_signalclass_in.signal');
signalclass_in.signal = signalclass_in.signal';
@@ -120,6 +120,9 @@ classdef EQ
% write to output
signalclass_out = signalclass_in;
noi = signalclass_out;
noi = signalclass_out - reference_signalclass_in;
end
function [yout,error_log] = process_(obj,data_in,ref_in)

33
Classes/04_DSP/Equalizer/FFE.m
View File

@@ -101,29 +101,6 @@ classdef FFE < handle
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;
@@ -153,16 +130,6 @@ classdef FFE < handle
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

14
Classes/04_DSP/Equalizer/FFE_DCremoval.m
View File

@@ -59,7 +59,7 @@ classdef FFE_DCremoval < handle
end
function [X] = process(obj, X, D)
function [X,Noi] = process(obj, X, D)
% actual processing of the signal (steps 1. - 3.)
% 1 normalize RMS
@@ -86,6 +86,8 @@ classdef FFE_DCremoval < handle
lbdesc = [num2str(obj.order),' tap FFE'];
X = X.logbookentry(lbdesc); % append to logbook
Noi = X;
Noi = X - D;
end
@@ -145,11 +147,21 @@ classdef FFE_DCremoval < handle
e_dc_buffer(1) = e_dc_est - obj.mu_dc * err(symbol);
e_dc_buffer = circshift(e_dc_buffer,1);
end
e_dc_save(symbol) = e_dc_est;
obj.error(epoch,symbol) = err(symbol) * err(symbol)'; % Instantaneous square error
end
end
if ~training
figure(1122)
hold on
scatter(1:numel(e_dc_save),e_dc_save,1,'.');
scatter(1:numel(y),y,1,'.');
end
end

35
Classes/04_DSP/Equalizer/FFE_FFDCAVG.m
View File

@@ -101,29 +101,6 @@ classdef FFE_FFDCAVG < handle
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
@@ -185,19 +162,7 @@ classdef FFE_FFDCAVG < handle
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

245
Classes/04_DSP/Equalizer/FFE_Kalman.m Normal file
View File

@@ -0,0 +1,245 @@
classdef FFE_Kalman < 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_Kalman(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,Noi] = 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
Noi = X;
Noi = X - D;
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)];
for epoch = 1 : epochs
symbol = 0;
% Initialization of Kalman filter variables
A = 1; % State transition matrix
H = 1; % Observation matrix
Q = 1e-4; % Process noise covariance
R = 1e-1; % Measurement noise covariance
P = 1; % Initial error covariance
mpi_est = 0; % Initial estimate for MPI noise
K = 0; % Kalman gain
subtract_mpi_est = 1;
for sample = 1 : obj.sps : N
symbol = symbol + 1;
% Get the current input sample and the equalizer output
U = x(obj.order + sample - 1 : -1 : sample);
if subtract_mpi_est || ~training
y(symbol,1) = (obj.e.' * U) - mpi_est .* 1 ; % Subtract MPI estimate
else
y(symbol,1) = obj.e.' * U;
end
% Decision and error calculation
if training
d_hat(symbol,1) = d(symbol);
else
[~, symbol_idx] = min(abs(y(symbol) - obj.constellation)); % Closest constellation point
d_hat(symbol,1) = obj.constellation(symbol_idx);
end
err(symbol) = y(symbol) - d_hat(symbol); % Instantaneous residual error
true_err(symbol) = y(symbol) - d(symbol);
% Kalman filter update to track the MPI noise
% Prediction step
P = A * P * A' + Q; % Update error covariance
K = P * H' / (H * P * H' + R); % Kalman gain
% Update step
mpi_est_new = mpi_est + K * (err(symbol) - H * mpi_est); % MPI noise estimation
alpha=0;
mpi_est = alpha * mpi_est + (1 - alpha) * mpi_est_new;
P = (1 - K * H) * P; % Update error covariance
% Subtract MPI noise from the signal
if subtract_mpi_est || ~training
y(symbol) = y(symbol);% - mpi_est;
else
end
% Equalizer weight update (LMS or NLMS)
if mio ~= 0
obj.e = obj.e - (mio * err(symbol) * U); % LMS weight update
else
normalizationfactor = (U.' * U);
obj.e = obj.e - err(symbol) * U / normalizationfactor; % NLMS weight update
end
% Store MPI estimate for visualization
mpi_estimates(symbol) = mpi_est;
Kgain(symbol) = K;
P_(symbol) = P;
H_(symbol) = H;
end %symbols
end %epoch
if ~training
if 1
% figure;
% subplot(2,2,1)
% hold on
% scatter(1:numel(y),y,1,'.');
% plot(1:numel(mpi_estimates), mpi_estimates);
% subplot(2,2,3)
% plot(1:numel(true_err), true_err);
% subplot(2,2,2)
% scatter(1:numel(y),y'-mpi_estimates,1,'.');
% xlabel('Sample Index');
% ylabel('MPI Noise Estimate');
% title('MPI Noise Estimation Over Time (Kalman Filter)');
% grid on;
figure(111);
subplot(3,1,1)
hold on
cla
scatter(1:numel(y),y,1,'.');
plot(1:numel(mpi_estimates), mpi_estimates,'DisplayName','mpi est');
ylim([-2,2]);
subplot(3,1,2)
cla
plot(1:numel(true_err), true_err,'DisplayName','true err');
ylim([-1,1]);
subplot(3,1,3)
plot(1:numel(true_err), true_err-mpi_estimates,'DisplayName',['diff rms: ',num2str(rms(true_err-mpi_estimates))]);
ylim([-1,1]);
legend
figure(333)
hold on
von = 5000;
bis = 15000;
plot(von:bis,true_err(von:bis),'DisplayName','Error')
plot(von:bis,movmean(true_err(von:bis), [30,30]),'DisplayName','Error')
plot(von:bis,mpi_estimates(von:bis),'DisplayName','Est','LineWidth',2);
legend
figure(222)
hold on
crr = xcorr(mpi_estimates,true_err,'normalized');
plot(crr);
end
% y = y-mpi_estimates';
end
end
end
end

197
Classes/04_DSP/Equalizer/FFE_Kalman_Feedback.m Normal file
View File

@@ -0,0 +1,197 @@
classdef FFE_Kalman_Feedback < 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_Kalman_Feedback(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,Noi] = 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);
% 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);
% 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
Noi = X;
Noi = X - D;
end
function [y, d_hat] = equalize(obj, x, d, mio, epochs, N, training)
arguments
obj
x
d
mio
epochs
N
training
end
% Initialize Kalman filter variables
A = 0.7; % State transition matrix
H = 1; % Observation matrix
Q = 1e-2; % Process noise covariance
R = 1e-3; % Measurement noise covariance
P = 1; % Initial error covariance
mpi_est = 0; % Initial estimate for MPI noise
K = 0; % Kalman gain
% Error buffers for parallelized structure
parallel_depth = 100; % Depending on your parallelization depth
err_buffer = zeros(parallel_depth, 1); % Store collected errors
mpi_est_buffer = zeros(parallel_depth, 1); % Store MPI estimates
x = [zeros(floor(obj.order/2), 1); x; zeros(obj.order, 1)];
for epoch = 1 : epochs
symbol = 0;
for sample = 1 : obj.sps : N
symbol = symbol + 1;
U = x(obj.order + sample - 1 : -1 : sample); % Input window for FFE
% Subtract the estimated MPI noise from the input signal before FFE
x_mpi_reduced = U;
y(symbol, 1) = obj.e.' * x_mpi_reduced; % Output of FFE
y(symbol, 1) = y(symbol, 1) - mpi_est;
% Decision and error calculation
if training
[~, symbol_idx] = min(abs(d(symbol) - obj.constellation)); % Closest constellation point
d_hat(symbol, 1) = d(symbol);
else
[~, symbol_idx] = min(abs(y(symbol) - obj.constellation)); % Closest constellation point
d_hat(symbol, 1) = obj.constellation(symbol_idx);
end
% Calculate the residual error
err(symbol) = y(symbol) - d_hat(symbol); % Instantaneous residual error
% Collect error for feedback after 'parallel_depth' samples
err_buffer(mod(symbol, parallel_depth) + 1) = err(symbol);
% Update Kalman filter based on the accumulated errors every parallel_depth samples
if mod(symbol, parallel_depth) == 0
% Compute the mean error over the last 'parallel_depth' samples
avg_error = mean(err_buffer);
% Prediction step (Kalman filter)
P = A * P * A' + Q; % Update the error covariance
K = P * H' / (H * P * H' + R); % Compute Kalman gain
% Update MPI noise estimate
mpi_est = mpi_est + K * (avg_error - H * mpi_est); % Update based on averaged error
P = (1 - K * H) * P; % Update error covariance
end
% Store MPI estimate for visualization or debugging
k_buffer(symbol) = K;
mpi_est_buffer(symbol) = mpi_est;
% FFE weight update using NLMS
if mio ~= 0
obj.e = obj.e - (mio * err(symbol) * U); % LMS weight update
else
normalizationfactor = (U.' * U);
obj.e = obj.e - err(symbol) * U / normalizationfactor; % NLMS weight update
end
end
end
if ~training
figure(111);
subplot(3,1,1)
hold on
cla
scatter(1:numel(y),y,1,'.');
plot(1:numel(mpi_est_buffer), mpi_est_buffer,'DisplayName','mpi est');
end
end
end
end

13
Classes/04_DSP/Equalizer/FFE_adaptive_decision.m
View File

@@ -55,7 +55,7 @@ classdef FFE_adaptive_decision < handle
end
function [X] = process(obj, X, D)
function [X,Noi] = process(obj, X, D)
% actual processing of the signal (steps 1. - 3.)
% 1 normalize RMS
@@ -80,11 +80,14 @@ classdef FFE_adaptive_decision < handle
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
Noi = X;
Noi = X - D;
end
function [y,d_hat] = equalize(obj,x,d,mio,epochs,N,training,showviz)
@@ -124,19 +127,21 @@ classdef FFE_adaptive_decision < handle
d_hat(symbol,1) = obj.constellation(symbol_idx);
end
y_buffer(symbol_idx,1) = y(symbol);
y_buffer(symbol_idx,:) = circshift(y_buffer(symbol_idx,:),1);
adap_constellation = mean(y_buffer,2,"omitnan");
delta_y_y(symbol) = y(symbol) - adap_constellation(symbol_idx);
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
obj.error(epoch,symbol) = err(symbol) * err(symbol)'; % Instantaneous square error