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

Question

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
x1=x+n;  % received noisy signal


%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?

Answer 1

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