clear
nr_data_bits=10000;% 0's and 1's, keep even number-Takes ~1 minute for a run of 1 million
nr_symbols=nr_data_bits/2;
b_data = (randn(1, nr_data_bits) > .5);%random 0's and 1's
b = (b_data);

d=zeros(1,length(b)/2);
%definition of the QPSK symbols using Gray coding.
for n=1:length(b)/2
    p=b(2*n);
    imp=b(2*n-1);
    if (imp==0)&&(p==0)
        d(n)=exp(j*pi/4);%45 degrees
    end
    if (imp==1)&&(p==0)
        d(n)=exp(j*3*pi/4);%135 degrees
    end
    if (imp==1)&&(p==1)
        d(n)=exp(j*5*pi/4);%225 degrees
    end
    if (imp==0)&&(p==1)
        d(n)=exp(j*7*pi/4);%315 degrees
    end
end
qpsk=d;

SNR=0:30;%change SNR values
BER1=[];
SNR1=[];
SER=[];
SER1=[];
sigma1=[];


for SNR=0:length(SNR);%loop over SNR-change SNR values (0,5,10 etc dB)
sigma = sqrt(10.0^(-SNR/10.0));
sigma=sigma/2;%Required a division by 2 to get close to exact solutions(Notes)-WHY?
x=randn(1,nr_symbols);
y=randn(1,nr_symbols);
ray=sqrt(0.5*(x.^2+y.^2));%variance=0.5-Tracks theoritical PDF closely
mpqpsk=qpsk.*ray;
mpsnqpsk=(real(mpqpsk)+sigma.*randn(size(mpqpsk))) +i.*(imag(mpqpsk)+sigma.*randn(size(mpqpsk)));


r=mpsnqpsk;%received signal plus noise and multipath
bhat=[real(r)<0;imag(r)<0];%detector
bhat=bhat(:)';
bhat1=bhat;%0's and 1's

ne=sum(b~=bhat1);%number of errors
BER=ne/nr_data_bits;
SER=ne/nr_symbols;%consider this to be Ps=log2(4)*Pb=2*Pb
SER1=[SER1 SER];
BER1=[BER1 BER];
SNR1=[SNR1 SNR];
end


figure(1);
plot(d,'o');%plot constellation without noise
axis([-2 2 -2 2]);
grid on;
xlabel('real'); ylabel('imag');
title('QPSK constellation');

figure(2);
semilogy(SNR1,BER1,'*',SNR1,SER1,'o');
grid on;
xlabel('SNR=Eb/No(dB)'); ylabel('BER/SER');
title('Simulation of BER/SER for QPSK with Gray coding( Rayleigh multipath and AWGN)');
legend('BER-simulated','SER-simulated');

figure(3)
plot(real(qpsk));
grid on;
axis([1 200 -2 2]);
title('QPSK symbols');
xlabel('symbols');ylabel('Amplitude');

figure(4)
plot(20*log10(abs(ray)));
grid on;
axis([1 200 -30 10]);
title('Rayleigh Fading Envelope(variance=0.5)');
xlabel('symbols');ylabel('Amplitude/RMS(dB)');