clear all;
close all;

%--------------- system parameters----------------
nt = 4;
nr = 4;
fft_len = 16; % FFt size
cp_len = 4; % cp length
L = 3; % channel order

F = 1/sqrt(fft_len)*exp(-j*2*pi/fft_len*[0:fft_len-1]'*[0:fft_len-1]); % FFT matrix
%--------------- simulation parameters-----------
SNR_loop = [0:5:30];
trial_loop = [50 50 50 50 50 50 50];
Pe = zeros(1,length(trial_loop));

%--------------- Monte Carlo --------------------
E = [];
for i = 1:length(SNR_loop)
SNR = SNR_loop(i);
NPW = 10^(-SNR/10);
trial = trial_loop(i);  
    for t = 1:trial %--------for loop for a specific SNR value
        %---------------- geneate source symbol----------------------
        data = randint(4,fft_len*2,2);   % 2 bits per symbol   
        data_MQAM = [];
        % QPSK modulation
        for m=1:fft_len
            data_MQAM = [data_MQAM 2*data(:,2*m-1) + data(:,2*m)];
        end
        % transmit data symbol
        S =[];
        for s = 1:4
            temp = modmap(data_MQAM(s,:),1,1,'qask',4);
            S = [S; 1/sqrt(2)*[temp(:,1)+j*temp(:,2)].'];   
        end    
            
        % equalvalent circulant matrix in time domaim
        %---------------- generate channel tap------------------------
        h  = 1/sqrt(2*nt*(L+1))*crandn(16,L+1);
        % equalvalent channel model in frequency domain
        H_freq = 1/sqrt(fft_len)*[fft(h.',16)].';
        H = zeros(nr,nt,nt*nr);
        S_r = [];
        for s = 1:fft_len  % tone by tone doing the flat fading MIMO process
            for n=1:nr
                H(n,:,s) = [H_freq((1+(n-1)*nt):n*nt,s).'];
            end
            S_r = [S_r H(:,:,s)*S(:,s) + (NPW/2)*crandn( nr, 1)]; 
        end
        
        %==================== receiver processing =========================
         
        %-------------- Per-tone MMSE VBLAST -------------------------------
        S_est = zeros(nr,fft_len);
        for s = 1:fft_len
            r_p = S_r(:,s);
            H_f = H(:,:,s);
            [S_est(:,s)] = VBLAST(H_f'*H_f, H_f , H_f'*r_p, r_p, nt, 1, NPW);
        end
        
        %-------------- error counting --------------------------------
        data = data.';
        data = data(:);
        S_est = S_est.';
        S_est = S_est(:);
        data_MQAM_est = demodmap(sqrt(2)*[real(S_est) imag(S_est) ],1,1,'qask',4);                      
        temp = dec2bin(data_MQAM_est,2);
        temp = temp.';
        temp = temp(:);
        data_est = [str2num(temp)].';  
        % caculate BER
        [numbit, error] = biterr(data.',data_est);
        Pe(i) = Pe(i) + error;

    end  % end trial
    Pe(i) = Pe(i)/trial;
    
end  % end SNR_loop
semilogy(SNR_loop,Pe)
xlabel('SNR dB')
ylabel('Pe')