352 lines
12 KiB
Matlab
352 lines
12 KiB
Matlab
classdef FFE_DCremoval_adaptive_mu < 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
|
|
e_tr
|
|
error
|
|
|
|
len_tr
|
|
mu_tr
|
|
epochs_tr
|
|
|
|
mu_dd
|
|
epochs_dd
|
|
|
|
mu_dc
|
|
dc_buffer_len
|
|
|
|
adaptive_mu_mode
|
|
|
|
ffe_buffer_len
|
|
|
|
smoothing_buffer_length
|
|
smoothing_buffer_update
|
|
|
|
constellation
|
|
|
|
decide
|
|
end
|
|
|
|
methods
|
|
function obj = FFE_DCremoval_adaptive_mu(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.mu_dc = 0.05;
|
|
options.dc_buffer_len = 1;
|
|
|
|
options.ffe_buffer_len = 1;
|
|
|
|
options.adaptive_mu_mode = 1;
|
|
|
|
options.smoothing_buffer_length = 0;
|
|
options.smoothing_buffer_update = 0;
|
|
options.decide = false;
|
|
|
|
end
|
|
|
|
assert(options.dc_buffer_len>0);
|
|
|
|
fn = fieldnames(options);
|
|
for n = 1:numel(fn)
|
|
obj.(fn{n}) = options.(fn{n});
|
|
end
|
|
|
|
obj.e = zeros(obj.order,1);
|
|
obj.error = 0;
|
|
|
|
obj.dc_buffer_len = floor(obj.dc_buffer_len);
|
|
|
|
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);
|
|
|
|
% if obj.smoothing_buffer_length > 0
|
|
% % Apply A1 filter smoothing
|
|
% % Calculate the moving sum with the window size N1
|
|
% moving_sum = movsum(X.signal, [obj.smoothing_buffer_length,0]);
|
|
%
|
|
% % Initialize the output smoothed signal
|
|
% X.signal = X.signal - (1 / obj.smoothing_buffer_length) * moving_sum;
|
|
% end
|
|
|
|
% Training Mode
|
|
training = 1;
|
|
obj.equalize(X.signal, D.signal,obj.mu_tr,obj.epochs_tr,obj.len_tr,training);
|
|
obj.e_tr = obj.e;
|
|
|
|
% Decision Directed Mode
|
|
N = X.length;
|
|
training = 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 - D;
|
|
|
|
end
|
|
|
|
function [y,d_hat] = equalize(obj, x, d, mu_lms, epochs, N, training)
|
|
% Equalize with adaptive DC-removal, VSS, and parallel-buffered DC updates
|
|
% Added: FFE gradient buffering in DD mode (error buffer) with update every obj.dc_buffer_len symbols
|
|
|
|
arguments
|
|
obj
|
|
x
|
|
d
|
|
mu_lms % LMS step-size (or 0 for NLMS)
|
|
epochs % number of training/DD epochs
|
|
N % number of samples to process
|
|
training % boolean flag: true->training mode, false->DD mode
|
|
end
|
|
|
|
if isempty(obj.e)
|
|
obj.e = zeros(obj.order,1);
|
|
end
|
|
|
|
% Zero-padding for filter memory
|
|
x = [zeros(floor(obj.order/2),1); x; zeros(obj.order,1)];
|
|
|
|
% Initialize storage
|
|
numSymbols = ceil(N/obj.sps);
|
|
y = zeros(numSymbols,1);
|
|
d_hat = zeros(numSymbols,1);
|
|
err = NaN(numSymbols,numel(obj.constellation));
|
|
e_dc_save= zeros(numSymbols,1);
|
|
|
|
% DC-adaptation parameters
|
|
P_err = 0; % running error power
|
|
alpha = 0.98; % forgetting factor for error power
|
|
err_prev = 0; % previous error sample for VSS correlation
|
|
gamma_dc = 1e-6; % meta step-size for DC VSS
|
|
mu_min = 1e-6; % lower bound for mu_dc
|
|
mu_max = 3e-1; % upper bound for mu_dc
|
|
|
|
% DC removal buffer
|
|
L = obj.dc_buffer_len; % buffer length
|
|
e_dc_buf = NaN(L,1);
|
|
e_dc_est = 0;
|
|
|
|
% FFE gradient buffer (DD mode only)
|
|
L_grad = obj.ffe_buffer_len; % buffer length
|
|
if ~training
|
|
% each column holds one past gradient of length obj.order
|
|
grad_buf = NaN(obj.order, L_grad);
|
|
end
|
|
|
|
smth_buffer = zeros(1, obj.smoothing_buffer_length);
|
|
smth_mean = 0;
|
|
% Main loop
|
|
for epoch = 1:epochs
|
|
s = 0;
|
|
for sample = 1:obj.sps:N
|
|
s = s + 1;
|
|
|
|
if obj.smoothing_buffer_length > 0
|
|
smth_buffer = circshift(smth_buffer,1,2);
|
|
smth_buffer(1) = x(sample);
|
|
if mod(s, obj.smoothing_buffer_update) == 0
|
|
smth_mean = mean(smth_buffer);
|
|
end
|
|
x(sample:sample+obj.sps-1) = x(sample:sample+obj.sps-1)-smth_mean;
|
|
end
|
|
|
|
U = x(obj.order+sample-1:-1:sample);
|
|
|
|
%-- 1) filter output with DC correction
|
|
y(s) = e_dc_est + obj.e.'*U;
|
|
|
|
%-- 2) decision
|
|
if training
|
|
[~, idx] = min(abs(d(s) - obj.constellation));
|
|
else
|
|
[~, idx] = min(abs(y(s) - obj.constellation));
|
|
end
|
|
d_hat(s) = obj.constellation(idx);
|
|
|
|
%-- 3) error
|
|
e_val = y(s) - d_hat(s);
|
|
if epoch == epochs
|
|
|
|
err(s,idx) = e_val;
|
|
true_err(s,idx) = y(s) - d(s);
|
|
end
|
|
|
|
%-- 4) tap-weight update: training immediate, DD buffered
|
|
if training
|
|
% immediate update (LMS or NLMS)
|
|
if mu_lms ~= 0
|
|
obj.e = obj.e - mu_lms * e_val * U;
|
|
else
|
|
normU = (U.'*U) + eps;
|
|
obj.e = obj.e - e_val * U / normU;
|
|
end
|
|
else
|
|
if 0
|
|
% buffer gradient
|
|
if mu_lms ~= 0
|
|
grad = e_val * U;
|
|
else
|
|
normU = (U.'*U) + eps;
|
|
grad = e_val * U / normU;
|
|
end
|
|
% shift and insert
|
|
grad_buf = circshift(grad_buf, 1, 2);
|
|
grad_buf(:,1) = grad;
|
|
% update once every L symbols
|
|
if mod(s, L_grad) == 0
|
|
avg_grad = mean(grad_buf, 2, 'omitnan');
|
|
if mu_lms ~= 0
|
|
obj.e = obj.e - mu_lms * avg_grad;
|
|
else
|
|
obj.e = obj.e - avg_grad;
|
|
end
|
|
end
|
|
end
|
|
end
|
|
|
|
|
|
%-- 5) DC adaptation
|
|
if obj.mu_dc ~= 0
|
|
|
|
if obj.adaptive_mu_mode
|
|
|
|
% VSS for mu_dc
|
|
delta_mu = gamma_dc * e_val * err_prev * (U.'*U);
|
|
obj.mu_dc = min(max(obj.mu_dc + delta_mu, mu_min), mu_max);
|
|
err_prev = e_val;
|
|
|
|
% DC buffer update & periodic estimate
|
|
P_err = alpha*P_err + (1-alpha)*e_val^2;
|
|
mu_dc_norm = obj.mu_dc / (P_err + eps);
|
|
|
|
else
|
|
|
|
% DC buffer update & periodic estimate
|
|
% P_err = alpha*P_err + (1-alpha)*e_val^2;
|
|
% mu_dc_norm = obj.mu_dc / (P_err + eps);
|
|
|
|
mu_dc_norm = obj.mu_dc;
|
|
|
|
end
|
|
|
|
e_dc_buf = circshift(e_dc_buf, 1);
|
|
e_dc_buf(1) = e_dc_est - mu_dc_norm * e_val;
|
|
|
|
if mod(s, L) == 0
|
|
e_dc_est = median(e_dc_buf, 'omitnan');
|
|
end
|
|
|
|
P_err_save(s) = P_err;
|
|
% Pcorr_save(s) = e_val * err_prev;
|
|
Ucorr_save(s) = (U.'*U);
|
|
mu_dc_save(s) = mu_dc_norm;
|
|
e_dc_save(s) = e_dc_est;
|
|
|
|
end
|
|
|
|
% store instantaneous squared error
|
|
obj.error(epoch, s) = e_val^2;
|
|
end
|
|
end
|
|
|
|
% Optional plotting in DD mode (uncomment if needed)
|
|
if 0%~training
|
|
|
|
constellation = unique(d);
|
|
lvlcol = cbrewer2('Paired', numel(constellation)*2);
|
|
lvlcol = lvlcol(2:2:end, :);
|
|
|
|
true_err(true_err==0) = NaN;
|
|
true_errmoverr = movsum(true_err, 4096, 'omitnan');
|
|
true_errmoverr = true_errmoverr./rms(true_errmoverr);
|
|
|
|
moverr = movsum(err, [100,100], 'omitnan');
|
|
moverr = moverr./rms(moverr);
|
|
|
|
figure(500); clf
|
|
hold on
|
|
% 1st subplot: true_errmoverr
|
|
% subplot(2,2,1); hold on
|
|
% for k = 1:4
|
|
% scatter(1:numSymbols, true_errmoverr(:,k), 1, lvlcol(k,:), '.');
|
|
% end
|
|
|
|
% scatter(1:numSymbols, Ucorr_save./rms(Ucorr_save), 1, lvlcol(1,:), '.','DisplayName','Ucorr_save');
|
|
% scatter(1:numSymbols, Pcorr_save./rms(Pcorr_save), 1, lvlcol(1,:), '.','DisplayName','P_corr');
|
|
% scatter(1:numSymbols, P_err_save, 1, lvlcol(1,:), '.','DisplayName','P_err');
|
|
% scatter(1:numSymbols, mu_dc_save, 1, lvlcol(2,:), '.','DisplayName','adapted value of $\mu_{DC}$');
|
|
% scatter(1:numSymbols, sum(moverr,2,'omitnan'), 1, lvlcol(1,:), '.','DisplayName','Mov Error $\hat{d}$ - x over all levels');
|
|
scatter(1:numSymbols, sum(e_dc_save,2,'omitnan'), 1, lvlcol(2,:), '.','DisplayName','Est. Error that is subtracted');
|
|
title('Moving Sum Error');
|
|
hold off
|
|
legend
|
|
|
|
% 2nd subplot: moverr
|
|
subplot(2,2,2); hold on
|
|
for k = 1:4
|
|
scatter(1:numSymbols, moverr(:,k), 1, lvlcol(k,:), '.');
|
|
end
|
|
title('Moving Sum Error');
|
|
hold off
|
|
legend
|
|
|
|
% 3rd subplot: err
|
|
subplot(2,2,3); hold on
|
|
for k = 1:4
|
|
scatter(1:numSymbols, err(:,k), 1, lvlcol(k,:), '.');
|
|
end
|
|
title('Error');
|
|
hold off
|
|
legend
|
|
|
|
% 4th subplot: err + obj.constellation'
|
|
subplot(2,2,4); hold on
|
|
for k = 1:4
|
|
scatter(1:numSymbols, err(:,k) + obj.constellation(k), 1, lvlcol(k,:), '.');
|
|
end
|
|
yline(obj.constellation, '--k');
|
|
title('Error + Constellation');
|
|
hold off
|
|
legend
|
|
|
|
sgtitle('Error Analysis Subplots');
|
|
|
|
end
|
|
end
|
|
|
|
|
|
end
|
|
end
|
|
|