clear;
close all;
symb_energyMin = 20;
symb_energyMax = 21;
symb_energy_d = 0.5;
packets_PerEnergystep = 20;
WGNPower = 1;
Transmitters = 2;
Receivers = 2;
M = input('qam?   ');
N=input('SM?   ');
packetSize =1000;%此處為符號數(shù)目
EnergiesdB = symb_energyMin:symb_energy_d:symb_energyMax;
Energies = 10.^(EnergiesdB./10);
for Energy_step = 1:size(Energies,2)
Eb =sqrt( Energies(Energy_step)/2);
    if N==1
            for bb=1:packets_PerEnergystep
            disp(['Eb = ',int2str(Energy_step),'; packet ',int2str(bb)]);
            H= randn (Receivers,Transmitters) + j*randn(Receivers,Transmitters);
            A = randint(Transmitters,packetSize,M); % Generate Data Stream A(rows =xmitters, cols = symbols)
            B = Eb*qammod(A,M); % 16-QAM Modulation
            noise = (WGNPower).*randn(size(H*B)) + j*(WGNPower).*randn(size(H*B));
            r = H*B + noise;
            G = pinv(H);
            QS=G*r/Eb;
            E=qamdemod(QS,M);
            SER(:,bb) = 1-sum(A==E,2)/size(A,2); % Measure BER per T-R pair.
            end

    else
            for bb=1:packets_PerEnergystep
            disp(['Eb = ',int2str(Energy_step),'; packet ',int2str(bb)]);
            A=randint(1,packetSize,M);
            B=zeros(2,packetSize);
            B(1,:)= Eb*qammod(A,M);
               for KK=0:499
                   B(2,2*KK+1)=B(1,2*KK+2);
                   B(1,2*KK+2)=-B(1,2*KK+2)';
                   B(2,2*KK+2)=B(1,2*KK+1)';
               end
            H= randn (Receivers,Transmitters) + j*randn(Receivers,Transmitters);
            noise = (WGNPower).*randn(size(H*B)) + j*(WGNPower).*randn(size(H*B));
             R = H*B + noise;
             E1=zeros(1,packetSize);
             for KK=0:499
                  E1(2*KK+1)=H(1,1)'* R(1,2*KK+1)+ H(1,2)* R(1,2*KK+2)'+H(2,1)'* R(2,2*KK+1)+H(2,2)'* R(2,2*KK+2)';
                  E1(2*KK+2)=H(1,2)'* R(1,2*KK+1)-H(1,1)* R(1,2*KK+2)'+H(2,2)'* R(2,2*KK+1)-H(2,1)* R(2,2*KK+2)';
             end
             E1=E1/(Eb*((abs(H(1,1)))^2 +(abs(H(1,2)))^2+(abs(H(2,1)))^2+(abs(H(2,2)))^2));
             E=qamdemod(E1,M);
             A-E
             
              SER(:,bb) = 1-sum(A==E,2)/size(A,2); 
            end
    end

p = mean(SER,2);
q = mean(p,1);
ERROR_RATE(Energy_step)=q;
end
figure;
semilogy(EnergiesdB,ERROR_RATE,'b-');
title (['\fontsize{12}\bfBLAST: \rm',int2str(M),' QAM']);
xlabel ('Eb/N_0'); ylabel ('SER');