Large MIMO system with ML detecor

Hi everybody please I need your help to implement Large MIMO system using V-blast technique. actually, I could to implement the code for ZF and MMSE for BPSK modulation. However, I can not achieve the ML detector when 4-QAM or 16-QAm are applied ( ZF and MMSE work FINE) ? i will be grateful if you can help me either by giving advice on my code or drive me to find the proper code in somewhere else.

*4-QAM code *

clear M=4; %Number of vectors in constellation %N=2048; %Number of samples (length of symbol)(Block size) N = 2048; % Block size total_errors=zeros(1,21); RandStream.setGlobalStream(RandStream('mt19937ar','seed', sum(100*clock))); %Random values for every session C=[-1+j,-1-j, 1-j, 1+j]; %4-QAM constellation vectors snr = [0:20]; % SNR values MT = 2; % Number of transmitters MR = 10; % Number of Receivers for i = 1:length(snr) for j=1:100 % p = rand(1,N)>0.5; % Generates random numbers % s = 2*p-1; % Converts 0 & 1 t0 -1 and 1 bk=randi([0 M-1],1,N);

        %Look-up table
        %Mod
        dk=C(bk+1);
sMod = kron(dk,ones(MR,1)); %
sMod = reshape(sMod,[MR,MT,N/MT]);
h = 1/sqrt(2)*[randn(MR,MT,N/MT) + j*randn(MR,MT,N/MT)]; % Rayleigh channel
n = 1/sqrt(2)*[randn(MR,N/MT) + j*randn(MR,N/MT)]; % White gaussian noise with 0dB variance
y = squeeze(sum(h.*sMod,2)) + 10^(-snr(i)/20)*n; % Noise addition
% forms weight vector for MMSE
hCof = zeros(2,2,N/MT) ;
hCof(1,1,:) = sum(h(:,2,:).*conj(h(:,2,:)),1) + 10^(-snr(i)/10);
hCof(2,2,:) = sum(h(:,1,:).*conj(h(:,1,:)),1) + 10^(-snr(i)/10);
hCof(2,1,:) = -sum(h(:,2,:).*conj(h(:,1,:)),1); % c term
hCof(1,2,:) = -sum(h(:,1,:).*conj(h(:,2,:)),1); % b term
hDen = ((hCof(1,1,:).*hCof(2,2,:)) - (hCof(1,2,:).*hCof(2,1,:)));
hDen = reshape(kron(reshape(hDen,1,N/MT),ones(2,2)),2,2,N/MT);
hInv = hCof./hDen;
hMod = reshape(conj(h),MR,N);
yMod = kron(y,ones(1,2));
yMod = sum(hMod.*yMod,1);
yMod = kron(reshape(yMod,2,N/MT),ones(1,2));
yHat = sum(reshape(hInv,2,N).*yMod,1);
pHat = real(yHat)>0; % Hard decision decoding

% nErr(i) = size(find([dk- pHat]),2); % Counts the # errors for k=1:N %% decision process [val,index]=min(abs(C-yHat(k))); bkest(k)=index-1; %best(i)=id-1;

            % BER computation

        end
        new_errors=biterr(bkest, bk);
        total_errors(i)=(total_errors(i)+new_errors);
    end   
     ber1(i)=total_errors(i)/(j*log2(M)*N);
end
% ber1 = nErr/N; % simulated ber

% BER(s)=total_errors(s)/(i*log2(M)*N); for i = 1:length(snr) for j=1:100 % p = rand(1,N)>0.5; % generating 0,1 with equal probability % s = 2*p-1; % BPSK modulation 0 -> -1; 1 -> 0 bk=randi([0 M-1],1,N);

        %Look-up table
        %Mod
        dk=C(bk+1);
sMod = kron(dk,ones(MR,1)); %
sMod = reshape(sMod,[MR,MT,N/MT]); % grouping in [MR,MT,N/MT ] matrix
h = 1/sqrt(2)*[randn(MR,MT,N/MT) + j*randn(MR,MT,N/MT)]; % Rayleigh channel
n = 1/sqrt(2)*[randn(MR,N/MT) + j*randn(MR,N/MT)]; % white gaussian noise, 0dB variance
y = squeeze(sum(h.*sMod,2)) + 10^(-snr(i)/20)*n;
hCof = zeros(2,2,N/MT) ;
hCof(1,1,:) = sum(h(:,2,:).*conj(h(:,2,:)),1);
hCof(2,2,:) = sum(h(:,1,:).*conj(h(:,1,:)),1);
hCof(2,1,:) = -sum(h(:,2,:).*conj(h(:,1,:)),1);
hCof(1,2,:) = -sum(h(:,1,:).*conj(h(:,2,:)),1);
hDen = ((hCof(1,1,:).*hCof(2,2,:)) - (hCof(1,2,:).*hCof(2,1,:)));
hDen = reshape(kron(reshape(hDen,1,N/MT),ones(2,2)),2,2,N/MT);
hInv = hCof./hDen;
hMod = reshape(conj(h),MR,N);
yMod = kron(y,ones(1,2));
yMod = sum(hMod.*yMod,1);
yMod = kron(reshape(yMod,2,N/MT),ones(1,2));
yHat = sum(reshape(hInv,2,N).*yMod,1);
pHat = real(yHat)>0;
for k=1:N
            %%decision process
            [val,index]=min(abs(C-yHat(k)));
            bkest(k)=index-1;
            %best(i)=id-1;

            % BER computation

        end
        new_errors=biterr(bkest, bk);
        total_errors(i)=(total_errors(i)+new_errors);
    end   
     ber2(i)=total_errors(i)/(j*log2(M)*N);
end
% ber1 = nErr/N; % simulated ber

% BER(s)=total_errors(s)/(i*log2(M)*N); for i = 1:length(snr) for j=1:100 bk=randi([0 M-1],1,N);

        %Look-up table
        %Mod
        dk=C(bk+1);
sMod = kron(dk,ones(MR,1));
sMod = reshape(sMod,[MR,MT,N/MT]);
h = 1/sqrt(2)*[randn(MR,MT,N/MT) + j*randn(MR,MT,N/MT)]; % Rayleigh channel
n = 1/sqrt(2)*[randn(MR,N/MT) + j*randn(MR,N/MT)]; % white gaussian noise, 0dB variance
y = squeeze(sum(h.*sMod,2)) + 10^(-snr(i)/20)*n;
% Maximum Likelihood Receiver
QAM_table = [-1+j,-1-j,  1-j, 1+j]/sqrt(10);
% if [s1 s2 ] = [+1,+1 ]
sHat1 = [1+j,-1+j];
sHat1 = repmat(sHat1,[1 ,N/2]);
sHat1Mod = kron(sHat1,ones(MR,1));
sHat1Mod = reshape(sHat1Mod,[MR,MT,N/MT]);
zHat1 = squeeze(sum(h.*sHat1Mod,2)) ;
J11 = sum(abs(y - zHat1),1);
% if [s1 s2 ] = [+1,-1 ]
sHat2 = [1+j -1-j];
sHat2 = repmat(sHat2,[1 ,N/2]);
sHat2Mod = kron(sHat2,ones(MR,1));
sHat2Mod = reshape(sHat2Mod,[MR,MT,N/MT]);
zHat2 = squeeze(sum(h.*sHat2Mod,2)) ;
J10 = sum(abs(y - zHat2),1);
% if [s1 s2 ] = [-1,+1 ]
sHat3 = [1+j 1-j];
sHat3 = repmat(sHat3,[1 ,N/2]);
sHat3Mod = kron(sHat3,ones(MR,1));
sHat3Mod = reshape(sHat3Mod,[MR,MT,N/MT]);
zHat3 = squeeze(sum(h.*sHat3Mod,2)) ;
J01 = sum(abs(y - zHat3),1);
% if [s1 s2 ] = [-1,-1 ]
sHat4 = [-1+j 1-j];
sHat4 = repmat(sHat4,[1 ,N/2]);
sHat4Mod = kron(sHat4,ones(MR,1));
sHat4Mod = reshape(sHat4Mod,[MR,MT,N/MT]);
zHat4 = squeeze(sum(h.*sHat4Mod,2)) ;
J00 = sum(abs(y - zHat4),1);
% if [s1 s2 ] = [-1,-1 ]
sHat5 = [-1+j -1-j];
sHat5 = repmat(sHat5,[1 ,N/2]);
sHat5Mod = kron(sHat5,ones(MR,1));
sHat5Mod = reshape(sHat5Mod,[MR,MT,N/MT]);
zHat5 = squeeze(sum(h.*sHat5Mod,2)) ;
J05 = sum(abs(y - zHat5),1);
% if [s1 s2 ] = [-1,-1 ]
sHat6 = [-1-j 1-j];
sHat6 = repmat(sHat6,[1 ,N/2]);
sHat6Mod = kron(sHat6,ones(MR,1));
sHat6Mod = reshape(sHat6Mod,[MR,MT,N/MT]);
zHat6 = squeeze(sum(h.*sHat6Mod,2)) ;
J06 = sum(abs(y - zHat6),1);

rVec = [J11;J10;J01;J00;J05;J06]; % Finds the minimum from the four possible combinations [u dd] = min(rVec,[],1); ref = [101 99; 101 -101; 101 -99; 99 -99 ; 99 -101; -101 -99]; % Bits mapping pHat = zeros(1,N); pHat(1:2:end) = ref(dd,1); pHat(2:2:end) = ref(dd,2); % [u ss] = min(bk,[],1); % mHat = zeros(1,N); % mHat = ref(u,1); % order=mHat.'; % mHat(2:2:end) = ref(ss,2); % nErr(i) = size(find([order-pHat]),2); gy=[-1+j,-1-j, 1-j, 1+j]; for k=1:N %% decision process [val,index]=min(abs(gy-pHat(k))); bkest(k)=index-1; %best(i)=id-1;

            % BER computation

        end
        new_errors=biterr(bkest, bk);
        total_errors(i)=(total_errors(i)+new_errors);
    end   
     ber3(i)=total_errors(i)/(j*log2(M)*N);
end
%  ber3 = nErr/N; % simulated ber
semilogy(snr,ber1,'r-',snr,ber2,'b-',snr,ber3,'g-'); % plots the ber for
legend('MMSE','ZF', 'ML');
xlabel('SNR (dB)');
ylabel('BER');
title('BER Performance of Different V-BLAST Detection Schemes over 2x10 MIMO Channels with Rayleigh Channel');
axis([0 20 10^-6 0.5])
grid on

16-QAM code

clear M=16; %Number of vectors in constellation %N=2048; %Number of samples (length of symbol)(Block size) N = 2048; % Block size total_errors=zeros(1,22); RandStream.setGlobalStream(RandStream('mt19937ar','seed', sum(100*clock))); %Random values for every session C=[-3+3*j, -3+j, -3-3*j, -3-j, -1+3*j, -1+j, -1-3*j, -1-j, 3+3*j, 3+j, 3-3*j, 3-j, 1+3*j, 1+j, 1-3*j, 1-j]; %16-QAM constellation vectors snr = [0:20]; % SNR values MT = 2; % Number of transmitters MR = 10; % Number of Receivers % for i = 1:length(snr) % for j=1:100 % % p = rand(1,N)>0.5; % Generates random numbers % % s = 2*p-1; % Converts 0 & 1 t0 -1 and 1 % bk=randi([0 M-1],1,N); % % %Look-up table % %Mod % dk=C(bk+1); % sMod = kron(dk,ones(MR,1)); % % sMod = reshape(sMod,[MR,MT,N/MT]); % h = 1/sqrt(2)*[randn(MR,MT,N/MT) + j*randn(MR,MT,N/MT)]; % Rayleigh channel % n = 1/sqrt(2)*[randn(MR,N/MT) + j*randn(MR,N/MT)]; % White gaussian noise with 0dB variance % y = squeeze(sum(h.*sMod,2)) + 10^(-snr(i)/20)*n; % Noise addition % % forms weight vector for MMSE % hCof = zeros(2,2,N/MT) ; % hCof(1,1,:) = sum(h(:,2,:).*conj(h(:,2,:)),1) + 10^(-snr(i)/10); % hCof(2,2,:) = sum(h(:,1,:).*conj(h(:,1,:)),1) + 10^(-snr(i)/10); % hCof(2,1,:) = -sum(h(:,2,:).*conj(h(:,1,:)),1); % c term % hCof(1,2,:) = -sum(h(:,1,:).*conj(h(:,2,:)),1); % b term % hDen = ((hCof(1,1,:).*hCof(2,2,:)) - (hCof(1,2,:).*hCof(2,1,:))); % hDen = reshape(kron(reshape(hDen,1,N/MT),ones(2,2)),2,2,N/MT); % hInv = hCof./hDen; % hMod = reshape(conj(h),MR,N); % yMod = kron(y,ones(1,2)); % yMod = sum(hMod.*yMod,1); % yMod = kron(reshape(yMod,2,N/MT),ones(1,2)); % yHat = sum(reshape(hInv,2,N).*yMod,1); % pHat = real(yHat)>0; % Hard decision decoding % % % % % nErr(i) = size(find([dk- pHat]),2); % Counts the # errors % for k=1:N % %% decision process % [val,index]=min(abs(C-yHat(k))); % bkest(k)=index-1; % %best(i)=id-1; % % % BER computation % % end % new_errors=biterr(bkest, bk); % total_errors(i)=(total_errors(i)+new_errors); % end % ber1(i)=total_errors(i)/(j*log2(M)*N); % end % % ber1 = nErr/N; % simulated ber % % % BER(s)=total_errors(s)/(i*log2(M)*N); % for i = 1:length(snr) % for j=1:100 % % p = rand(1,N)>0.5; % generating 0,1 with equal probability % % s = 2*p-1; % BPSK modulation 0 -> -1; 1 -> 0 % bk=randi([0 M-1],1,N); % % %Look-up table % %Mod % dk=C(bk+1); % sMod = kron(dk,ones(MR,1)); % % sMod = reshape(sMod,[MR,MT,N/MT]); % grouping in [MR,MT,N/MT ] matrix % h = 1/sqrt(2)*[randn(MR,MT,N/MT) + j*randn(MR,MT,N/MT)]; % Rayleigh channel % n = 1/sqrt(2)*[randn(MR,N/MT) + j*randn(MR,N/MT)]; % white gaussian noise, 0dB variance % y = squeeze(sum(h.*sMod,2)) + 10^(-snr(i)/20)*n; % hCof = zeros(2,2,N/MT) ; % hCof(1,1,:) = sum(h(:,2,:).*conj(h(:,2,:)),1); % hCof(2,2,:) = sum(h(:,1,:).*conj(h(:,1,:)),1); % hCof(2,1,:) = -sum(h(:,2,:).*conj(h(:,1,:)),1); % hCof(1,2,:) = -sum(h(:,1,:).*conj(h(:,2,:)),1); % hDen = ((hCof(1,1,:).*hCof(2,2,:)) - (hCof(1,2,:).*hCof(2,1,:))); % hDen = reshape(kron(reshape(hDen,1,N/MT),ones(2,2)),2,2,N/MT); % hInv = hCof./hDen; % hMod = reshape(conj(h),MR,N); % yMod = kron(y,ones(1,2)); % yMod = sum(hMod.*yMod,1); % yMod = kron(reshape(yMod,2,N/MT),ones(1,2)); % yHat = sum(reshape(hInv,2,N).*yMod,1); % pHat = real(yHat)>0; % for k=1:N % %% decision process % [val,index]=min(abs(C-yHat(k))); % bkest(k)=index-1; % %best(i)=id-1; % % % BER computation % % end % new_errors=biterr(bkest, bk); % total_errors(i)=(total_errors(i)+new_errors); % end % ber2(i)=total_errors(i)/(j*log2(M)*N); % end % ber1 = nErr/N; % simulated ber % ber1 = nErr/N; % simulated ber

% BER(s)=total_errors(s)/(i*log2(M)*N); for i = 1:length(snr) for j=1:100 bk=randi([0 M-1],1,N);

        %Look-up table
        %Mod
        dk=C(bk+1);
sMod = kron(dk,ones(MR,1));
sMod = reshape(sMod,[MR,MT,N/MT]);
h = 1/sqrt(2)*[randn(MR,MT,N/MT) + j*randn(MR,MT,N/MT)]; % Rayleigh channel
n = 1/sqrt(2)*[randn(MR,N/MT) + j*randn(MR,N/MT)]; % white gaussian noise, 0dB variance
y = squeeze(sum(h.*sMod,2)) + 10^(-snr(i)/20)*n;
% Maximum Likelihood Receiver
% if [s1 s2 ] = [+1,+1 ]
% sHat1 = [1+j,-1+j];
% sHat1 = repmat(sHat1,[1 ,N/2]);
% sHat1Mod = kron(sHat1,ones(MR,1));
% sHat1Mod = reshape(sHat1Mod,[MR,MT,N/MT]);
% zHat1 = squeeze(sum(h.*sHat1Mod,2)) ;
% J11 = sum(abs(y - zHat1),1);
% % if [s1 s2 ] = [+1,-1 ]
% sHat2 = [1+j -1-j];
% sHat2 = repmat(sHat2,[1 ,N/2]);
% sHat2Mod = kron(sHat2,ones(MR,1));
% sHat2Mod = reshape(sHat2Mod,[MR,MT,N/MT]);
% zHat2 = squeeze(sum(h.*sHat2Mod,2)) ;
% J10 = sum(abs(y - zHat2),1);
% % if [s1 s2 ] = [-1,+1 ]
% sHat3 = [1+j 1-j];
% sHat3 = repmat(sHat3,[1 ,N/2]);
% sHat3Mod = kron(sHat3,ones(MR,1));
% sHat3Mod = reshape(sHat3Mod,[MR,MT,N/MT]);
% zHat3 = squeeze(sum(h.*sHat3Mod,2)) ;
% J01 = sum(abs(y - zHat3),1);
% % if [s1 s2 ] = [-1,-1 ]
% sHat4 = [-1+j 1-j];
% sHat4 = repmat(sHat4,[1 ,N/2]);
% sHat4Mod = kron(sHat4,ones(MR,1));
% sHat4Mod = reshape(sHat4Mod,[MR,MT,N/MT]);
% zHat4 = squeeze(sum(h.*sHat4Mod,2)) ;
% J00 = sum(abs(y - zHat4),1);
% % if [s1 s2 ] = [-1,-1 ]
% sHat5 = [-1+j -1-j];
% sHat5 = repmat(sHat5,[1 ,N/2]);
% sHat5Mod = kron(sHat5,ones(MR,1));
% sHat5Mod = reshape(sHat5Mod,[MR,MT,N/MT]);
% zHat5 = squeeze(sum(h.*sHat5Mod,2)) ;
% J05 = sum(abs(y - zHat5),1);
% % if [s1 s2 ] = [-1,-1 ]
% sHat6 = [-1-j 1-j];
% sHat6 = repmat(sHat6,[1 ,N/2]);
% sHat6Mod = kron(sHat6,ones(MR,1));
% sHat6Mod = reshape(sHat6Mod,[MR,MT,N/MT]);
% zHat6 = squeeze(sum(h.*sHat6Mod,2)) ;
% J06 = sum(abs(y - zHat6),1);
% 
% rVec = [J11;J10;J01;J00;J05;J06]; % Finds the minimum from the four possible combinations
% [u dd] = min(rVec,[],1);
% ref = [101 99; 101 -101; 101 -99; 99 -99 ; 99 -101; -101 -99]; % Bits mapping
% pHat = zeros(1,N);
% pHat(1:2:end) = ref(dd,1);
% pHat(2:2:end) = ref(dd,2);
% % nErr(i) = size(find([d- pHat]),2);
% gy=[-1+j,-1-j,  1-j, 1+j];
QAM_table = [-3-3j, -3-j, -3+3j, -3+j, -1-3j, -1-j, -1+3j, -1+j,
3-3j, 3-j, 3+3j, 3+j, 1-3j, 1-j, 1+3j, 1+j]/sqrt(10);
metric = 100000;
tt=sum(h.*sMod,2);
sHat1Mod =  reshape(tt,MR,N/MT);
for l = 1:16
    x1_tmp=QAM_table(l);
x_tmp(1) = QAM_table(l); Esti_y1 = y - sHat1Mod.*x1_tmp;
for m = 1:16
     x2_tmp=QAM_table(m);
x_tmp(2) = QAM_table(m);
Esti_y2 = Esti_y1 - sHat1Mod.*x2_tmp;
for n = 1:16
     x3_tmp=QAM_table(n);
    x_tmp(3) = QAM_table(n);
Esti_y3 = Esti_y2 - sHat1Mod.*x3_tmp;
for o = 1:16
     x4_tmp=QAM_table(o);
x_tmp(4) = QAM_table(o);
Esti_y4 = Esti_y3 - sHat1Mod.*x4_tmp;
metric_tmp = sqrt(Esti_y4'*Esti_y4);
final_dis=reshape (Esti_y4,5,N);
if metric_tmp<metric
X_hat = x_tmp; metric = metric_tmp;
end
end
end
end
end
for k=1:N
            %%decision process
            [val,index]=min(abs(C-Esti_y4(k)));
            bkest(k)=index-1;
            %best(i)=id-1;

            % BER computation

        end
        new_errors=biterr(bkest, bk);
        total_errors(i)=(total_errors(i)+new_errors);
    end   
     ber3(i)=total_errors(i)/(j*log2(M)*N);
     % % nErr(i) = size(find([dk- pHat]),2); % Counts the # errors

end % ber3 = nErr/N; % simulated ber semilogy(snr,ber3,'g-'); % plots the ber for legend('MMSE','ZF', 'ML'); xlabel('SNR (dB)'); ylabel('BER'); title('BER Performance of Different V-BLAST Detection Schemes over 2x10 MIMO Channels with Rayleigh Channel'); axis([0 20 10^-6 0.5]) grid on

many THanks