# Matlab: Help in applying equalizer and proper calculation of symbol error rate for QAM modulation in Rayleigh channel

I am unable to understand why the Constant Modulus ALgorithm is not functioning properly for QAM modulation in Rayleigh multipath fading. The same code worked for BPSK in AWGN (previous Question asked Matlab: What should be the BER performance for BPSK using Constant Modulus Algorithm equalizer.)

The problem is in the equalization step, the variable yy is containing NaN values and so the weights are not getting calculated. What should I do to make it work? Thank you.

UPDATE 1 : I have updated the code with the new algorithm explained in the answer but still the problem persists which is the algorithm does not converge. I am unsure if I am calculating the symbol error rate properly or not. Would need help in whether to calculate SER or BER

UPDATE 2 : I am confused regarding the analytical expression for symbol error rate and bit error rate. When I used the analytical formula given in Book : Andrea Goldsmith and the function berawgn() in Matlab-- both give different results.

 % Script for computing the BER for QAM modulation in Rayleigh
% channel

clear all
%N = 1e5; % number of bits or symbols
% Script for computing the BER for QAM modulation in Rayleigh
% channel

clear all
%N = 1e5; % number of bits or symbols
Eb_N0_dB = [0:15]; % multiple Eb/N0 values

nTaps = 3;
mu = 0.001;

j = sqrt(-1);
ht = 1/sqrt(2)*[randn(1,nTaps) + j*randn(1,nTaps)];

L = length(ht);

M=64;
N = 10^5;
% Convert Eb/No values to channel SNR
% modulation
k = sqrt(1/((2/3)*(M-1))); % normalizing factor
m = [1:sqrt(M)/2]; % alphabets
alphaMqam = [-(2*m-1) 2*m-1];

j = sqrt(-1);
Eb_N0_dB = Eb_N0_dB + 10*log10(k);

for ii = 1:length(Eb_N0_dB)

uncoded_bits  = round(rand(1,N));

data=ceil(M.*rand(N,1))-1;
tx=((randi(8,1,N)*2-9)+j*(randi(8,1,N)*2-9))/sqrt(42);  %64 QAM

ht = 1/sqrt(2)*[randn(1,nTaps) + j*randn(1,nTaps)];

%rx = tx.*ray;
rx = filter(ht,1,tx);  %multi path Rayleigh fading
% Noise variance
N0 = 1/10^(Eb_N0_dB(ii)/10);

rx = rx + sqrt(N0/2)*(randn(1,length(tx))+1i*randn(1,length(tx)));

%CMA
Le = 5; %Equalizer length
e = zeros(length(rx),1); % error
w = ones(Le,1); % equalizer coefficients
%w(Le)=1; % actual filter taps are flipud(w)!
yd = zeros(length(rx),1);
xx = zeros(Le, 1);
R2 = 1.380953; %64 QAM ; 1.32 for 16 QAM
yy = w(:,i).'*xx;
yd(i)= yy;
e(i) = (yy)*(R2-abs(yy)^2);  %MODIFIED ALGORITHM ERROR FUNCTION
w(:,i+1) = w(:,i) + mu*e(i)*conj(yy); %MODIFIED
%e(i) = abs(yy)^2 - R2;
%w(:,i+1) = w(:,i) - mu * xx * e(i) * yy.';
end
sb = sb';
% M-QAM demodulator at the Receiver

% demodulation
y_re = real(sb)/k;
y_im = imag(sb)/k;
% rounding to the nearest alphabet
% 0 to 2 --> 1
% 2 to 4 --> 3
% 4 to 6 --> 5 etc
ipHat_re = 2*floor(y_re/2)+1;
ipHat_re(find(ipHat_re>max(alphaMqam))) = max(alphaMqam);
ipHat_re(find(ipHat_re<min(alphaMqam))) = min(alphaMqam);

% rounding to the nearest alphabet
% 0 to 2 --> 1
% 2 to 4 --> 3
% 4 to 6 --> 5 etc
ipHat_im = 2*floor(y_im/2)+1;
ipHat_im(find(ipHat_im>max(alphaMqam))) = max(alphaMqam);
ipHat_im(find(ipHat_im<min(alphaMqam))) = min(alphaMqam);

ipHat = ipHat_re + j*ipHat_im;

BER(ii) = size(find(data - ipHat'),2);

end

AWGNtheoryBer = 0.5*erfc(sqrt(10.^(Eb_N0_dB/10))); % theoretical ber

EbN0Lin = 10.^(Eb_N0_dB/10);

Pe_sim=[];
for i = 1:length(Eb_N0_dB),
Pe_sim=[Pe_sim BER(i)];
end

% plotting the simulated data

% comparing with theoretical results

if(M==4)
a=1;
else
a=4/log2(M);
end
b=3*log2(M)/(M-1);
RaytheoryBer=[];
for i = 1:length(Eb_N0_dB),
RaytheoryBer= [RaytheoryBer 0.5*a*(1-sqrt(0.5*b*EbN0Lin(i)/(1+0.5*b*EbN0Lin(i))))];
end                      % Ref: Wireless Communication, A. Goldsmith

figure;
semilogy(Eb_N0_dB,Pe_sim,'r-','Linewidth',2);
hold on;
semilogy(Eb_N0_dB, AWGNtheoryBer,'bs-');
hold on;
semilogy(Eb_N0_dB, RaytheoryBer,'g-');
hold on;

semilogy(Eb_N0_dB, SERTheoryMatlab,'k-');
axis([Eb_N0_dB(1) Eb_N0_dB(end) 10^-5 0.5])
grid on
legend('sim-CMA','Theoretical-BER-AWGN channel','Raytheory-Goldsmith','Raytheory-Matlab');
xlabel('Eb/No, dB');
ylabel('Bit Error Rate');
title('Bit error probability curve for QAM in ISI with Rayleigh Channel');


Let's take a look at the cost function of CMA:

$J(\mathbf{w}) = E[|y_{k}^{2} - 1|^{2}]$,

onde $y_{k} = \mathbf{w}^{H} \mathbf{x}_{k}$ and $\mathbf{x}_{k}$ is the filter input.

It says that the update function will converge to coefficients that approximate the output of the filter to something whose squared value is $1$. So, it is appropriate to use CMA when you have in the input a BPSK signal, because the symbols in BPSK are ${-1,1}$, and its squared values are $1$.

But, it will not work if your equalizer input are not signals whose squared values could be 1.

I don't know the problem with your code, but in the foirst moment, CMA will not work for QAM modulation.

I have never worked with Multi Modulus Algorithm (MMA), but I think that this would be an appropriate algorithm to apply insted of CMA.

UPDATE:

According to the paper you posted in the comments, he is apllying a variation of the CMA, that have the following cost function (equation (5) in the paper):

$J(\mathbf{w}) = E[(|x_{k}|^{2} - \Delta_{2})^{2}]$,

where $\Delta_{2} = E[|x(k)|^{4}]/E[|x(k)|^{2}]$. That leads to the following error function (equation (6) in the paper):

$e(k) = x(k)(\Delta_{2} - |x_{k}|^{2})$

And then, the update function becomes (equation (7) in the paper):

$w(i+1) = w(i) + \mu e(k) v^{*}(k - i)$

So, in your code you are implementing the CMA that has the first cost function I wrote, and not the second one. You need to implement the second one in your code.

Read this paper again, carefully, and you will be able to implement by yourself.

• Thank you for pointing out the reason. I found a paper radioeng.cz/fulltexts/2011/11_03_683_691.pdf where on pg 12, Fig 12 gives the simulation for 16QAM. So, in order to make CMA work for 16QAM could you please help in what I should do? I am not getting the curve as given in the paper. The plot will help me to justify the need for MMA. Thank you for your effort and time once again. – Ria George Nov 12 '15 at 16:24
• As I said, the cost function as I wrote above, cannot be applied to high order QAM modulation. You need to make changes in the cost function in order to use this kind of modulation. The author of the paper that you linked did it, but didn't show how. And I don't know, also. You will need to make a research about this topic. – JohnMarvin Nov 12 '15 at 19:37
• Actually, He showed, I am sorry. I will update my answer. – JohnMarvin Nov 13 '15 at 0:20
• Facing a lot of difficulty in implementing the Author's version – Ria George Nov 13 '15 at 4:08
• I went through the paper and the equations that you have plotted. I have Questions as how to calculate the modulus of the signal, Δ2. Should I be doing Xf=fft(x); Xf=fftshift(Xf); modulus_R = abs(Xf)^4 / abs(Xf)^2 ?here x is the filter output – Ria George Nov 13 '15 at 4:28