Fractionally spaced equalizer using LMS algorithm

Question

I am doing a project about the fractionally spaced equalizers on matlab. The goal is to mitigate ISI on the signal using the LMS algorithm to update the tap weights of the FSE. Below i have uploaded the whole code. I have done most of the project but have some problems on the final part. When i try to write the LMS algorithm part it seems to always have some problems (below i've put 2 versions of it, one is commented %%). I don't know what might be the problem. Since i'm using a BPSK modulated signal the final result should be a graph with points scattered around -1 and 1 forming 2 clouds arond these 2 values. Can you help me understand where does my error lies and how to fix it? The help will be very much appreciated.

Below i have explained all the code to make it easier for everyone to follow:

Channel impulse response creation. (RRC channel impulse response).

% Root of raised cosine filter parameters
alpha = 0.35;    % Rolloff factor
L     = 24;      % Truncated impulse length (-L/2:L/2)
Nc    = 12;      % Number of samples per symbol 

% Tx filter design
istrt = floor(L/2);
n     = -istrt:1/Nc:istrt;     %(-L/2:1/Nc:L/2)
pt    = truncRRC(n, alpha,0);

% Plot Tx/Rx impulse response
figure, stem(n, pt), title('RRC filter impulse response')

Insertion of ISI in the channel impulse response. Since i have to have ISI in my channel impulse response (otherwise i wouldn't need equalization) i have created an identical copy of "pt" called "pt_delay" that is just shifted in time. After that i have used this method to insert ISI "pt_delay = 0.7pt+0.3pt_delay;"

pt_delay = truncRRC(n, alpha,1/3);
pt_delay = 0.7*pt+0.3*pt_delay;
figure, stem(n, pt_delay), title('RRC filter impulse response with delay')

Here the BPSK modulated signal is created.

% M-PSK consellation
M      = 2;                           % Number of PSK levels
mb     = log2(M);                     % Number of bits
psklev = exp(1j*2*pi*(0:M-1)/M);      % M-PSK levels
nSym   = 1000;                        % Number of symbols

symtx       = randi([0,M-1], nSym, 1);  % Generate random M-ary symbols
[symtxg, ~] = togray(symtx, mb);        % Gray mapping
ci          = psklev(symtxg+1);         % Symbol selection

After that i have created my packet that consists of a Barker preamble, training sequence and the BPSK symbols. Upsample the signal and then filter it.

pairs = [1 2; 3 4; 5 6; 7 8; 9 10; 11 12; 13 14;
         15 16; 17 18; 19 20; 21 22; 23 24; 25 26;
         27 28 ;29 30; 31 32; 33 34; 35 36; 37 38;
         39 40; 41 42; 43 44; 45 46; 47 48; 49 50;
         51 52; 53 54; 55 56; 57 58; 59 60; 61 62; 63 64];

% Generate Barker sequence
barker = ones(1, 16);
for i = 1:size(pairs, 1)
    barker(pairs(i, 1)) = -1;
    barker(pairs(i, 2)) = 1;
end
barker_len = size(barker);

% Generate training sequence
trainSymtx  = randi([0, M-1], 1, nSym);   % Random M-ary symbols

ReadySignal = [barker,trainSymtx, ci];

% Transmission filtering: upsample and filter
txSig_up = upsample(ReadySignal, Nc);
txSig    = filter(pt_delay, 1, txSig_up);

Noise addition

SNR    = 30;                          % SNR in dB
SNRlin = 10^(SNR/10);                 % SNR linear scale
swn    = sqrt(0.5*Nc/(SNRlin*mb));    % Noise variance
rxSig  = txSig + swn*(randn(size(txSig)) + 1j*randn(size(txSig)));

figure, scatter(real(rxSig), imag(rxSig)), title('Received constellation')

% Receiver filtering
rxSig2 = filter(pt, 1, rxSig);

After adding noise and filtering it's time for cross correlation to find where does my "useful" signal starts.

preamble_up = upsample(barker, Nc);
[r,lags]    = xcorr(rxSig2,preamble_up);  

[~, idx]           = max(abs(r));
start_preamble_idx = lags(idx);    % start_preamble_idx is the index of the symbol that precedes the first of the preamble

I have downsampled to pass from 12 samples per symbol to only 2 samples per symbol

% downsampling
rxSig_ds = downsample(rxSig2(start_preamble_idx+1:end), Nc/2);

figure, scatter(real(rxSig_ds), imag(rxSig_ds)), title('Received constellation downsampled')

Here i have divided my received packet into its 3 parts

% ML detection (minimum distance)
dist_vec           = abs(psklev.' - rxSig_ds) .^ 2;
[~, sym_idx]       = min(dist_vec);
det_symg           = sym_idx - 1;                   % Gray-coded detected symbol
[detected_syms, ~] = fromgray(det_symg, mb);        % Detected symbol
selected_syms      = psklev(detected_syms + 1);     % Symbol selection

% select only the preamble
preamble_syms = selected_syms(1:2*barker_len(2));

% select only the training sequence
trainseq_syms = selected_syms(2*barker_len(2)+1:2*barker_len(2)+2*nSym);

% select only the payload
payload_syms = selected_syms(2*barker_len(2)+2*nSym+1:end);

figure, scatter(real(payload_syms), imag(payload_syms)), title('Received payload')
ylim([-0.5 0.5])

Now we can check if the vector 'preamble_syms' contains the initial preamble 'barker'.

preamble_syms_real = real(preamble_syms);
preamble_syms_real = downsample(preamble_syms_real, 2);

if preamble_syms_real == barker
    disp('The vectors are equal.');
else
    disp('The vectors are not equal.');
end

Equalizer start

% FSE parameters
tapSpacing = 1/Nc;  % Fractional tap spacing
numTaps    = L*Nc;  % Number of taps

% LMS algorithm parameters
stepSize = 0.01;   

% Initialize tap weights
tapWeights = zeros(1, numTaps);

train_syms_real    = real(trainseq_syms);
error              = zeros(1, numTaps);
train_syms_real_ds = downsample(train_syms_real, 2);

This is where LMS algorithm is applied and where my problems start. As you can see i have tried 2 different approaches but nothing seems to work.

% % Perform equalization for each symbol in the training signal
% for i = 1:2*length(trainSymtx)
%     
%     % Extract the current training symbol
%     currentSymbol = train_syms_real(i);
%     
%     % Convolve the tap weights with the received training seq using the fractional tap positions
%     convOutput = conv(tapWeights, currentSymbol);
%     
%     % Calculate the error
%     error = trainSymtx(i) - convOutput;
%     
%     % Update the tap weights using the LMS update rule
%     tapWeights = tapWeights + stepSize * conj(currentSymbol) * error;
% 
% end
for i = 1:nSym

    convOutput(i) = tapWeights * train_syms_real_ds';
    %convOutput(i) = conv(tapWeights, train_syms_real_ds);
    error(i) = trainSymtx(i)-convOutput(i);

    for m = 1:nSym

        tapWeights(m) = tapWeights(m) + stepSize * error(i) * train_syms_real_ds(m);
    end

end

% % Use the trained tap weights to equalize the received signal
% equalizedSignal = conv(tapWeights, payload_syms);

% %Downsample to have 1 sample per symbol
% equalizedSignal_ds = downsample(equalizedSignal,2);


figure, scatter(real(equalizedSignal_ds), imag(equalizedSignal_ds)), title('Constellation after equalization')

Necessary functions

function pRRC=truncRRC(x,a,tau)
%truncRRC: Truncated raised root cosine impulse
%Use:   pRRC=truncRRC(x,a,tau)
%       x: normalised time
%       a: rolloff, within [0,1]
%       tau: input delay
%       truncRRC(x,a)= (1-a) sinc( (1-a) x) +
%       a( sinc(ax + 1/4) cos(pi (x+1/4)) + sinc(ax - 1/4) cos(pi (x-1/4))

x = x - tau; % delayed time
pRRC = (1-a)*sinc(x*(1-a));
pRRC = pRRC + a* sinc(a*x + 0.25).*cos(pi * (x + 0.25));
pRRC = pRRC + a* sinc(a*x - 0.25).*cos(pi * (x - 0.25));
end

function [gd, gb] = togray(in, len)

if ischar(in)
in = bin2dec(in);
end

gd = bitxor(in, bitshift(in,-1));
gb = dec2bin(gd, len);

end

function [de, bi] = fromgray(in, len)

if ischar(in)
in = bin2dec(in);
end

in = in(:);
n_syms = length(in);

b = zeros(n_syms,len);

% get MSB of in
b(:,1) = bitget(in, repmat(len,[n_syms,1]));

for sh = 2 : len
    in_loc = bitget(in, repmat(len-sh+1,[n_syms,1]));
    b(:,sh) = bitxor(b(:,sh-1), in_loc);
end

bi = num2str(b);
de = bin2dec(bi);
end

Answer 1

There's a lotta content and a lotta detail. I just skipped over to the LMS adaptive filter.

for i = 1:nSym

    convOutput(i) = tapWeights * train_syms_real_ds';
    %convOutput(i) = conv(tapWeights, train_syms_real_ds);
    error(i) = trainSymtx(i)-convOutput(i);

    for m = 1:nSym

        tapWeights(m) = tapWeights(m) + stepSize * error(i) * train_syms_real_ds(m);
    end

end

Now it seems to me that the number of iterations of i is the length of the signal and, for a real-time process is effectively $\infty$. But the number of iterations of m is the number of taps of your FIR filter, much smaller than for 'i'.

Then the only other suggestion is to really investigate the value of stepSize. You might also want to consider the Normalized LMS adaptive filter. Unnormalized, your adaptation speed or rate is proportional to stepSize and the power of the input and the adaption might become unstable for loud signals. Normalized LMS will make the adaptation rate (to the tap coefficients) the same whether the general signal environment is loud or quiet.

If I recall, the LMS filter is ($x[n]$ input, $d[n]$ desired, $y[n]$ output, $e[n]$ error):

$$\begin{align} y[n] &= \sum\limits_{i=0}^{I_\text{order}} h[i,n] \, x[n-i] \\ \\ e[n] &= d[n] - y[n] \\ \\ h[i,n+1] &= h[i,n] \ - \ g_\text{step} \, e[n] \, x[n-i] \qquad 0 \le i \le I_\text{order}\\ \end{align}$$

The Normalized LMS will replace the last line with something that will reduce $g_\text{step}$ by something reciprocal to the current sliding mean of $|x[n]|^2$. I don't remember a good way to do that, probably square and sliding sum is the simplest cheap way, or maybe a first-order IIR (leaky integrator) for the sliding sum. I have to think about it, I haven't done this for 3 decades.

Maybe this: $$\begin{align} y[n] &= \sum\limits_{i=0}^{I_\text{order}} h[i,n] \, x[n-i] \\ \\ e[n] &= d[n] - y[n] \\ \\ h[i,n+1] &= h[i,n] \ - \ g_\text{step} \frac{ e[n] \, x[n-i] }{\overline{|x|^2}[n]+\epsilon} \qquad 0 \le i \le I_\text{order} \\ \\ \overline{|x|^2}[n+1] &= p \cdot \overline{|x|^2}[n] + (1-p) |x[n]|^2\\ \end{align}$$