function z=vblast_est(alg,modulation,corr,alpha,system)

%**************************************************************************
% This program implements 1x2 and 2x2 Vblast based on preamble type of
% channel estimation using MMSE estimation, on a frame of 130 symbols.
% z=vblast_est(alg,modulation,corr,alpha,system) where
% alg -> algorithm for vblast receiver.Choose from 'ZF' (for zero-forcing),
% 'MM'(MMSE estimation)
% modulation -> 'BPSK,'QPSK','16QAM','64QAM' 
% corr -> 1 for correlation at the receiver(2x2 correlation only implemented),
% and 0 for no correlation
% alpha -> correlation coefficient value from 0 -> 1 and 0 in case corr = 0;
% system -> 12 for a 1x2 system and 22 for a 2x2 system
% EXAMPLE: vblast_est('ZF','QPSK',1,0.5,22)plots ZF curve for QPSK with 0.5  
% correlation coefficient at the receiver and channel estimation for a 2x2
% system
%**************************************************************************

% identify bits for each type of modulation
switch modulation   
    case 'BPSK'
        BITS=1;
    case 'QPSK'
        BITS=2;
    case '16QAM'
        BITS=4;
    case '64QAM'
        BITS=6;
end

%define frame
K=130;
% define SNR range
EbNo=[0:2:30]; 

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

%set initial count
idx=1; 

%define numbers of antennas
if system==12
    M=2;%rx antennas
    N=1;%tx antennas
elseif system==22
    M=2;%rx antennas
    N=2;%tx antennas
end

%number of monte-carlo runs
Num=10;

%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 & system==22
                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(K,N,BITS)>0;   
            temp1=[];
                for i=1:K   
                    d1=tx_modulate(tx_bits(i,:),modulation);
                    temp1=[temp1; d1];
                end
                d=temp1;
%insert pilots
        if system==12
            pilots=[0 0 0 0 0 0 0 0 0 0 1 -1 -1  1 -1 -1  1 -1  1  1];
            d=[pilots.'; temp1];
        elseif system==22
            pilots=[0 0 0 0 0 0 0 0 0 0 1 -1 -1  1 -1 -1  1 -1  1  1;1 -1 -1 -1  1 -1 -1  1 -1  1 0 0 0 0 0 0 0 0 0 0];
            d=[pilots.'; temp1];
        end

            
%AWGN noise            
         	AWGN_noise = sqrt(sigma)*(randn(150, M)+j*randn(150, 2));
%receiver signal vector added to AWGN noise            
            r = ((H_save*d.')/sqrt(N)).' + AWGN_noise;
           
            pilots=pilots.';
            tr_symbols=r(1:20,:);
            R=r(21:end,:).';
        if system==12
            h11=mean(tr_symbols(:,1).*conj(pilots(:,1)));
            h12=mean(tr_symbols(:,2).*conj(pilots(:,1)));
            est_coefs=[h11;h12];
        elseif system==22
            h11=mean(tr_symbols(:,1).*conj(pilots(:,1)));
            h12=mean(tr_symbols(:,1).*conj(pilots(:,2)));
            h21=mean(tr_symbols(:,2).*conj(pilots(:,1)));
            h22=mean(tr_symbols(:,2).*conj(pilots(:,2)));
            est_coefs=[[h11;h21]  [h12; h22]];
        end
            H_save=est_coefs;
for i=1:K
r=    R(:,i);
X=temp1.';
X=X(:,i);
H=H_save;
            

%Zero-forcing algorithm            
            if alg=='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=='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
            end
%count errors    
            if BITS==1 %for BPSK modulation
                errors(iter) =  sum((sign(real(a))~=sign(real(X))));
            else %for QPSK,16QAM and 64QAM modulations
                errors(iter) =  sum((sign(real(a))~=sign(real(X))) | sign(imag(a))~=sign(imag(X)));
            end
        end %end of iteration loop Loop
    end
    BER(idx)=sum(errors)/(Num) ; % Calculate BER after completion of 'Num' runs
    SER(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');