function z=vblast(Num,alg,modulation,corr,alpha)

%**************************************************************************
% This program implements vblast using perfect channel estimation for a 2x2
% system.
% z=vblast(Num,alg,modulation,corr,alpha) where
% Num-> number of runs
% alg -> algorithm.Choose from 'ZF' (for zero-forcing),'MM'(MMSE
% estimation),'ML'(maximum-likelihood decoding)
% modulation -> 'BPSK,'QPSK','16QAM','64QAM' (Note: for 'ML' algorithm only
% BPSK and QPSK exist and that too only for 2x2 configuration, as otherwise 
% run time is too high.
% corr -> 1 for correlation at the receiver, 0 for no correlation
% alpha -> correlation coefficient value from 0 -> 1 and 0 in case corr =0;
% EXAMPLE: vblast(1000,'ZF','QPSK',1,0.5)plots ZF curve for QPSK with 0.5  
% correlation coefficient at the receiver using 1000 Monte Carlo runs
%**************************************************************************

clc;
close all;
% identify bits for each type of modulation
SEED = 931316785;
rand('seed',SEED);

Num=1000;
alg=['ZF' 'MM' 'ML'];
modulation='QPSK';
corr=0;
alpha=0;
switch modulation   
    case 'BPSK'
        BITS=1;
    case 'QPSK'
        BITS=2;
    case '16QAM'
        BITS=4;
    case '64QAM'
        BITS=6;
end
for class=1:3
    
% define SNR range
EbNo=[0:2:30]; 

%define plotting axis
SNR_axis=[]; 
BER_axis=[]; 

%set initial count
idx=1; 

%define numbers of antennas
M=2;%rx antennas
N=2;%tx antennas


%parameters for wait bar
h = waitbar(0,'Please wait...');
wb=6.25;

%clear BER register
BER=[];

%commence SNR loop
for SNR=EbNo 
    errors=0;
%define standard deviation, sigma,  for noise     
    sigma=0.5/(sqrt(N)*10^(SNR/10)); 
        for iter=1:Num %commence iteration loop
%define a random Rayleigh channel            
            H=(randn(M,N)+j*randn(M,N))/sqrt(2);
%exercise correlation option, if any
            if corr
                R=chol([1 alpha;alpha 1]);%cholesky factor of correlation matrix, i.e. 'R' square root of correlation matrix
                H=R*H; % R^0.5 *H -> correlation at receiver
            end
            H_save=H;%assign H to unalterable value

% modulated input data
            tx_bits=randn(N,BITS)>0;   
            temp1=[];
                for i=1:N   
                    d1=tx_modulate(tx_bits(i,:),modulation);
                    temp1=[temp1; d1];
                end
            d=temp1;
%AWGN noise            
         	AWGN_noise = sqrt(sigma)*(randn(M, N)+j*randn(M, N));
%receiver signal vector added to AWGN noise            
            r = (H_save*d)/sqrt(N) + sqrt(sigma)*(randn(M, 1)+j*randn(M, 1));
%Zero-forcing algorithm            
            if alg(class)=='ZF'
 %initialization
                G=pinv(H);
                [gk k0]=min(sum(abs(G).^2,2));
      
                    for m=1:N     % This FOR loop determines the ordering sequence k1 and determines the 'a' matrix. 
                                  %This is just one run,i.e. one for each H matrix.
                 	    		  %The 'a' matrix is automatically sorted as [a1 a2...aM]
                        k1(m)=k0;
                        w(m,:)=G(k1(m),:);
                        y=w(m,:)*r;
                        a(k1(m),1)=Q(y,modulation);
                        r = r - a(k1(m)) * H_save(:, k1(m));   
                        H(:,k0)=zeros(M,1);
                        G=pinv(H);
                    for t=1:m
                        G(k1(t),:)=inf;
                    end
                    [gk k0]=min(sum(abs(G).^2,2));
                end
%MMSE algorithm                
            elseif alg(class)=='MM'
  %initialisation
                    G=inv(H'*H+N/(10^(0.1*SNR))*eye(N))*H';
                    [gk k0]=min(sum(abs(G).^2,2));
        
      
                    for m=1:N     % This FOR loop determines the ordering sequence k1 and determines the 'a' matrix.  
                                  % This is just one run,i.e.one for each H matrix.The 'a' matrix is automatically sorted 
                                  % as [a1 a2...aM]
         
                        k1(m)=k0;
                        w(m,:)=G(k1(m),:);
                        y=w(m,:)*r;
                        a(k1(m),1)=Q(y,modulation);
                        r = r - a(k1(m)) * H_save(:, k1(m));   
                        H(:,k0)=zeros(M,1);
                        G=inv(H'*H+N/(10^(0.1*SNR))*eye(N))*H';
                        for t=1:m
                            G(k1(t),:)=inf;
                        end
                        [gk k0]=min(sum(abs(G).^2,2));
                    end
%ML algorithm                    
            elseif alg(class)=='ML'
                    p=form_ref_matrix(BITS); %create file containing 2x2 sets of constellation symbols as reference
                    temp2=[];
                    temp3=H*p/sqrt(N);
                    if BITS==1
                        temp4=4; %square(number of symbols in constellation (2 for BPSK)), since this is a 2x2 system
                    elseif BITS==2
                        temp4=16;%square(number of symbols in constellation (4 for QPSK)), since this is a 2x2 system
                    end
                    for i=1:temp4
                        temp2(:,i)=abs(r-temp3(:,i)).^2;
                    end
                    w=sum(temp2);
                    [y1 x1]=min(w);
                    a=[p(1,x1); p(2,x1)];
                    temp2=[];
                    w=[];
            end
%count errors    
            if BITS==1 %for BPSK modulation
                errors(iter) =  sum((sign(real(a))~=sign(real(d))));
            else %for QPSK,16QAM and 64QAM modulations
                errors(iter) =  sum((sign(real(a))~=sign(real(d))) | sign(imag(a))~=sign(imag(d)));
            end
        end %end of iteration loop Loop
    BER(o,idx)=sum(errors)/(Num) ; % Calculate BER after completion of 'Num' runs
    SER(o,idx)=BER(idx)*BITS; %calculate symbol error rate
    idx=idx + 1; %increment count
    waitbar(wb/100);
    wb=wb+6.25;%increment wait bar
end %end of SNR loop
close(h);%terminate wait bar

SNR_axis=EbNo;
BER_axis=[BER_axis BER];
SER_axis=SER;

%plot BER
semilogy(SNR_axis,BER_axis,'b-*');
xlabel('SNR [dB]');
ylabel('BER/SER');
title('BER/SER Plots');

hold;

%plot SER
semilogy(SNR_axis,SER_axis,'b-o');
axis([0 30 1e-6 1]);
grid on;
legend('BER','SER');

hold off;
if alg=='MM'
    alg='MMSE'
end
str=['VBLAST System-' '2 x 2 ' alg ' Algorithm with ' modulation ' Modulation'];
set(gcf,'NumberTitle','off');
set(gcf,'Name',str);
grid on