Matlab: Error when working with higher order QAM signal - Matrix dimension must agree

Question

This problem seems to be trivial but I am left scratching my head when trying to resolve it. I am trying to apply Fractionally spaced equalizer with constant modulus technique for 64 QAM constellation. The program works for QPSK or 4 QAM but when I apply it to 64QAM, it throws error :

Error using  / 
Matrix dimensions must agree.

Error in Working_FSE_CMA_64QAM (line 68)
sb1=sb/(fh(temp));  % scale the output

I don not have the Communications toolbox so have generated 64QAM symbols using the answer given in my previous question Generate 16 QAM signal

Can somebody please help in making the code work? Thank you.

% Blind channel estimation/equalization
% adpative CMA method in Fractional space

T=1000;    % total number of data
dB=25;     % SNR in dB value

%%%%%%%%% Simulate the Received noisy Signal  %%%%%%%%%%%
N=5; % smoothing length N+1
Lh=5;  % channel length = Lh+1
Ap=4;  % number of subchannels or receive antennas

h=randn(Ap,Lh+1)+sqrt(-1)*randn(Ap,Lh+1);   % channel (complex)
for i=1:Ap, h(i,:)=h(i,:)/norm(h(i,:));    end        % normalize
s = (randi(8,1,T)*2-5)+j*(randi(8,1,T)*2-5);  %64 QAM
%s=round(rand(1,T))*2-1;  % QPSK or 4 QAM symbol sequence
%s=s+sqrt(-1)*(round(rand(1,T))*2-1);

% generate received noisy signal
x=zeros(Ap,T);    % matrix to store samples from Ap antennas
SNR=zeros(1,Ap);
for i=1:Ap
    x(i,:)=filter(h(i,:),1,s);
    vn=randn(1,T)+sqrt(-1)*randn(1,T);   % AWGN noise (complex)
    vn=vn/norm(vn)*10^(-dB/20)*norm(x(i,:));  % adjust noise power
    SNR(i)=20*log10(norm(x(i,:))/norm(vn));   % Check SNR of the received samples
    x(i,:)=x(i,:)+vn;                        % received signal
end
SNR=SNR    % display and check SNR

%%%%%%%%%%%%% adaptive equalizer estimation via CMA
Lp=T-N;   %% remove several first samples to avoid 0 or negative subscript
X=zeros((N+1)*Ap,Lp);  % sample vectors (each column is a sample vector)
for i=1:Lp
    for j=1:Ap
        X((j-1)*(N+1)+1:j*(N+1),i)=x(j, i+N:-1:i).';
    end
end

e=zeros(1,Lp);  % used to save instant error
f=zeros((N+1)*Ap,1); f(N*Ap/2)=1;    % initial condition
%R2=2;                  % constant modulas of QPSK symbols
R2 = 1.380953; %For 64 QAM http://www.google.com/patents/US7433400
mu=0.001;      % parameter to adjust convergence and steady error
for i=1:Lp
   e(i)=abs(f'*X(:,i))^2-R2;                  % instant error
   f=f-mu*2*e(i)*X(:,i)*X(:,i)'*f;     % update equalizer 
   f(N*Ap/2)=1;
%   i_e=[i/10000 abs(e(i))]             % output information 
end   

%sb=f'*X;   % estimate symbols (perform equalization)
sb = filter(f, 1, X);
% calculate SER
H=zeros((N+1)*Ap,N+Lh+1);  temp=0;
for j=1:Ap
    for i=1:N+1, temp=temp+1; H(temp,i:i+Lh)=h(j,:); end  % channel matrix
end

fh=f'*H;    % composite channel+equalizer response should be delta-like 
temp=find(abs(fh)==max(abs(fh)));   % find the max of the composite response

sb1=sb/(fh(temp));  % scale the output
sb1=sign(real(sb1))+sqrt(-1)*sign(imag(sb1));  % perform symbol detection
start=N+1-temp;  % general expression for the beginning matching point
sb2=sb1(10:length(sb1))-s(start+10:start+length(sb1));  % find error symbols
SER=length(find(sb2~=0))/length(sb2)   % calculate SER

if 1
    subplot(221), 
    plot(s,'o');   % show the pattern of transmitted symbols
    grid,title('Transmitted symbols');  xlabel('Real'),ylabel('Image')
    axis([-2 2 -2 2])

    subplot(222),
    plot(x,'o');  % show the pattern of received samples
    grid, title('Received samples');  xlabel('Real'), ylabel('Image')

    subplot(223),
    plot(sb,'o');   % show the pattern of the equalized symbols
    grid, title('Equalized symbols'), xlabel('Real'), ylabel('Image')

    subplot(224),
    plot(abs(e));   % show the convergence
    grid, title('Convergence'), xlabel('n'), ylabel('Error e(n)')
end

Answer 1

The main problem is that the iterative algorithm does not converge (or more specifically diverges). As a results the values in f contains NaN. The result of

temp=find(abs(fh)==max(abs(fh)))

is then an empty vector which thus does not match the expected scalar size on the line

sb1=sb/(fh(temp));  % scale the output

To fix the problem, you may note that the value of mu and R2 you use are based on a signal constellation with unit variance. You can generate such a 64-QAM constellation with:

s = ((randi(8,1,T)*2-9)+j*(randi(8,1,T)*2-9))/sqrt(42);  %64 QAM

Also, a few lines redefine the complex constant j=sqrt(-1) to use as an index variable. You should avoid using j as index. Also, clearing all variables at the beginning of your sript (with clear all) can help start you execution with a consistent clean state.