# MATLAB : Proper estimation of weights and how to calculate MSE for QPSK signal for Constant Modulus Algorithm

I had asked conceptual Question about Constant Modulus Algorithm. I am implementing the simple steps of the algorithm, as I do not have the cma() built in module.

I am considering a Finite Impulse Response (FIR) system whose true coefficients are $h = \begin{bmatrix}1 &0.45 &-0.2\end{bmatrix}$. The algorithm is presented briefly :

$u_n$ : Output of the FIR filter that is driven by Input signal, $s$ is a white Gaussian signal that drives an FIR process by the command

s=round(rand(1,N))*2-1;
u = filter(h,1,s);


We receive a noisy signal which is the output of the system corrupted by AWGN. Let the noise corrupted signal be

$$x_n = h^T s_n + \eta_n$$ Constructing an equalizer

$$y_n = w^T x_n$$

The cost function to be minimized by gradient descent is

$$J(w) = E\left[\left(\big\lvert y_n\big\rvert^2 -1\right)^2\right]$$

The weight update equation is given by

$$w_{n+1} = w_n - 2\mu e_n y_n^T x_n$$

where error,

$$e_n = \left(\big\lvert y_n\big\rvert^2 -1\right)$$

Operator $T$ is the transpose and assuming real signals without any imaginary and complex parts.

Below is the code. The algorithm does not give proper result in terms of estimates of the weights.

Problem 1: The Graph for the Mean Square Error (in DB) for the estimated weights returned by the algorithm vs. Signal to Noise Ratio is giving an opposite trend i.e. instead of MSE decreasing with increasing SNR, I am getting the opposite when I re-run the program!! Below is the image of what I mean. The first figure is correct, but with the same code, I ran it again and got the second figure. Why is this & how can I prevent it?

I am not sure if this is due to the initialization of the variable $L$, filter order = number of delays = 2. Somethings are not straight for me which are for the FIR filter of the form : $u(t) = e(t) + 0.45e(t-1) - 0.2e(t-2)$, what is the filter order, smoothing length and the number of weights unknown to be estimated? Can somebody please help so that it works?

   clear all
clc

N = 256;
h = [1    0.45   -0.2];

R2 = 2;

mu =  1.0000e-009;
noisedB  =0;

L=2; % smoothing length L+1
ChL=1;  % length of the channel= ChL+1
EqD=round((L+ChL)/2);  %  channel equalization delay

i=sqrt(-1);

Ch=[1 0.45 -0.2]; %Channel
%Ch=[0.8+i*0.1 .9-i*0.2]; %complex channel
Ch=Ch/norm(Ch);% normalize

skip =1
for l=1:6
i          = 1;

TxS=round(rand(1,N))*2-1;  % QPSK symbols are transmitted symbols
%TxS=TxS+sqrt(-1)*(round(rand(1,N))*2-1);
x=filter(Ch,1,TxS); %channel distortion

n=randn(1,N);   % additive white gaussian noise

n=n/norm(n)*10^(-noisedB/20)*norm(x);  % scale noise power

%estimation using CMA
K=N-L;   %% Discard initial samples for avoiding 0's and negative
X=zeros(L+1,K);  %each vector
for j=1:K
X(:,j)=x1(j+L:-1:j).';  %y_n = w^T x_n
end

e=zeros(1,K);
w=zeros(L+1,1);
w(EqD)=1;    % initial condition
while i<=K

e(i)=abs(w.'*X(:,i))^2-R2;                  % initial error
w=w-mu*2*e(i)*X(:,i)*X(:,i)'*w;  % update equalizer co-efficients
cma_mse_h(l,i)    = sum((w'-h).^2)/3;
est_w(i,:) = w;
w(EqD)= 1;

i = i+1;

end
noisedB = noisedB + 5;
end

for ii = 1:6

Error(ii) = 10*log10(mean(cma_mse_h(ii,:)));
end

plot([0:5:25], Error(1:6));
grid on;
xlabel ('SNR(dB)')
ylabel('MSE_h')


Second Problem : The true FIR channel coefficients are not imaginary and have only real parts. But the estimated weights will have both real & imaginary if I work with real and complex representation. How can I properly calculate the weights & its MSE for this case?

• Some things to try: experiment with your step size, try using normalized LMS instead, and try increasing the number of taps in your equalizer. – Jason R May 21 '15 at 12:25
• @JasonR: I played around with the step size but not the delay. I am not sure whether for blind CMA we assume the knowledge of the order of the filter. But, yet the same misconvergence issue remains. Another thing is that I am not quite clear what the variable ChL, EqD stand for : I have modified the CMA part of the code found in link mathworks.com/matlabcentral/fileexchange/… Can you please help me with these initial problems? Thank you for your suggestions. – SKM May 21 '15 at 15:47

Your code reveals many misconception about what the CMA is supposed to achieve:

• your step size mu is much too small; note, however, that the optimal step size can only be found through experiment.
• the variabe noisedB appears to be the desired SNR of the received signal. An SNR of $0\,\text{dB}$ as specified by you is very poor (the noise is as strong as the signal), and in such noisy conditions the CMA cannot perform properly.
• you have a channel with $3$ taps, and your equalizer length is also only $3$ taps; this is much too short to achieve any reasonable equalization.
• you define the MSE as the error between the channel impulse response and the equalizer response. This is wrong, because the equalizer is not supposed to converge to the channel impulse response! The equalizer must equalize the channel, i.e. the concatenation of the channel and the equalizer should give a pure delay.

There are a few other small errors in your code, but instead of correcting them all I show you a very simple and short code example. The figure below the code shows the error curve vs symbol index, and the total impulse response of the channel and the equalizer. As you can see, it approximates a delta impulse (i.e. no distortion). The SNR is $30\,\text{dB}$ for this example. The error becomes much smaller if there is no noise. For worse SNRs, the algorithm becomes unstable.


N = 20000;          % # symbols
h = [1,.45,-.2];    % channel impulse response
h = h/norm(h);
Le = 20;            % equalizer length
mu = .001;          % step size
snr = 30;           % snr in dB

s0 = round( rand(N,1) )*2 - 1;     % BPSK signal
s = filter(h,1,s0);                % filtered signal

% add Gaussian noise at desired snr
n = randn(N,1);
vs = var(s);
vn = vs*10^(-snr/10);
n = sqrt(vn)*n;
r = s + n;          % received signal

e = zeros(N,1);     % error
w = zeros(Le,1);    % equalizer coefficients
w(Le)=1;            % actual filter taps are flipud(w)!

yd = zeros(N,1);

for i = 1:N-Le,
x = r(i:Le+i-1);
y = w'*x;
yd(i)=y;
e(i) = y^2 - 1;
w = w - mu * e(i) * y * x;
end

np = 100;           % # sybmols to plot (last np will be plotted); np < N!

subplot(3,1,1), plot(e.*e), title('error')
subplot(3,1,2), stem(conv(flipud(w),h)), title('equalized channel impulse response')
subplot(3,1,3), plot(1:np,s0(N-np+1:N),1:np,yd(N-np+1-Le+1:N-Le+1))
title('transmitted and equalized signal'), legend('transmitted','equalized'), axis([0,np,-1.5,2])



• Thank you for your detailed clarification and the code. Could you shed some guidelines on these as these are the objectives that I wanted to achieve hence I had used the MSE in my code ? (1) The thing is I need to compare the performance of CMA in blind system identification and equalization by using the Cramer Rao bound of the coefficients for blind system identification with the MSE of the channel impulse (estimated coefficients). – SKM May 24 '15 at 6:23
• Theoretically, MSE of the estimated coefficients $w$ should asymptotically reach the CRLB with increasing SNR. That is the reason I was applying the MSE formula in my code. How can I obtain the MSE of the weights/ coefficients $w$? – SKM May 24 '15 at 6:23
• @SKM: (1) This is a totally different question. Your original question is about equalizing a channel, which is different from system identification. So for this point, you should formulate a new question, which clearly states what you need to know. (2) You should compare the equalized signal to the original transmitted signal before the channel, because the equalizer tries to compensate for the channel. (3) The filter order is indeed 2 (order = length-1). – Matt L. May 24 '15 at 7:16
• @SKM: That's right, $y$ is the equalized signal which should approximate the originally transmitted signal. – Matt L. May 24 '15 at 10:10
• @SKM: The second plot is the total impulse response of the channel and the equalizer. If it is a unit impulse then the equalizer compensates perfectly for the channel distortion. As you can see from the first plot (error over time), after 2000 symbols the equalizer has not yet fully converged, so the total impulse response is not yet a unit impulse. As for the plot of the original and equalized symbols, I've added the corresponding code and plot to my answer. – Matt L. May 25 '15 at 17:49