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);

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));

        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)]; 

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

        % Send over Gaussian Link to the receiver
        rx = rx + sqrt(N0/2)*(randn(1,length(tx))+1i*randn(1,length(tx)));

        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);
        received_sig = rx';
        xx = zeros(Le, 1);
        R2 = 1.380953; %64 QAM ; 1.32 for 16 QAM
        for i = 1:length(received_sig)
            xx = [received_sig(i); xx(1:end-1)];
            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.';
        sb = filter(w(:,length(received_sig)+1), 1,  received_sig); % estimate symbols (perform equalization)
        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);


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

EbN0Lin = 10.^(Eb_N0_dB/10);

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

% plotting the simulated data

% comparing with theoretical results

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

    [BERRaytheorymatlab SERTheoryMatlab] = berfading(Eb_N0_dB,'qam',64,1);

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');

1 Answer 1


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.


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.

  • $\begingroup$ 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. $\endgroup$
    – Ria George
    Commented Nov 12, 2015 at 16:24
  • $\begingroup$ 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. $\endgroup$
    – JohnMarvin
    Commented Nov 12, 2015 at 19:37
  • $\begingroup$ Actually, He showed, I am sorry. I will update my answer. $\endgroup$
    – JohnMarvin
    Commented Nov 13, 2015 at 0:20
  • $\begingroup$ Facing a lot of difficulty in implementing the Author's version $\endgroup$
    – Ria George
    Commented Nov 13, 2015 at 4:08
  • $\begingroup$ 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 $\endgroup$
    – Ria George
    Commented Nov 13, 2015 at 4:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.