% MPSK
%
% K. Bell
% 11/22/99

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Simulation Parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

M = 8;
K = log2(M);
switch M
case 2                    % BPSK symbols
   sym_map=[1;-1];
   Es = 10.^([[-17:3:-2] [-1:1:13]]/10); % Es
case 4                    % QPSK symbols
   sym_map=[1;j;-1;-j];
   Es = 10.^([[-17:3:1] [2:1:16]]/10); % Es
case 8                    % 8PSK symbols
   sym_map=[1;(1+j)/sqrt(2);j;(-1+j)/sqrt(2);-1;(-1-j)/sqrt(2);-j;(1-j)/sqrt(2)];
   Es = 10.^([[-17:3:7] [8:1:22]]/10); % Es
end

Eb = Es/K;
na = length(Es);
Ns = 100000;                % Number of symbols
No = 2;                   % noise unit variance

Es_No = Es/No;
Eb_No = Eb/No;

bits = round(rand(K,Ns));           % KxNs matrix of random 0,1 bits
symbols = graymapPSK(bits);         % gray code map to symbols
BER = zeros(1,na);
SER = zeros(1,na);
Pseint = zeros(1,na);

dphi = 0.01*pi/M;
phi = [-pi/M+dphi/2:dphi:pi/M];
nphi = length(phi);
dv = 0.01;

for m=1:na
   m
   r = sqrt(Es(m))*symbols + sqrt(No/2)*randn(2,Ns);         % observations
   cr = r(1,:)+j*r(2,:);
   sd = zeros(2,Ns);
   for n=1:Ns
      [ee ind]=min(abs(sym_map-cr(n)));                    % map to closest symbol
      sd(:,n)= [real(sym_map(ind));imag(sym_map(ind))];    % symbol decision
   end
   bd = grayunmapPSK(sd,M);
   errors = abs(bd-bits);                % KxNs matrix of bit errors
   symb_err = sign(sum(errors,1));       % one symbol error per column
   SER(m) = sum(symb_err)/Ns;
   BER(m) = sum(sum(errors))/(K*Ns);
   p1 = zeros(1,nphi);
   for q = 1:nphi
      mv = sqrt(2*Es_No(m))*cos(phi(q));
      v = [max(dv/2,mv-5+dv/2):dv:mv+5];
      pv = exp(-0.5*((v-mv).^2))/sqrt(2*pi);
      p1(q) = sum(v.*pv)*dv; 
   end
   ep = exp(-Es_No(m)*sin(phi).^2).*p1/sqrt(2*pi);
   Pseint(m) = 1-sum(ep)*dphi;

end

% erfc*(x) = 0.5*erfc(x/sqrt(2))

switch M
case 2
   
   Pse = 0.5*erfc(sqrt(2*Eb_No)/sqrt(2));
   Pbe = Pse;

case 4
   Pbe = 0.5*erfc(sqrt(2*Eb_No)/sqrt(2));
   Pse = 1-(1-Pbe).^2;

case 8
   Pse = Pseint;
   Pbe = Pse/K;
end

Pslb = 0.5*erfc(sqrt(2*Es_No)*sin(pi/M)/sqrt(2));
Psub = 2*Pslb;


figure(1);
subplot(1,2,1)
   semilogy(10*log10(Es_No),SER,'oc')
   hold on
   semilogy(10*log10(Es_No),Pse,'g')
   semilogy(10*log10(Es_No),Pseint,'b')
   semilogy(10*log10(Es_No),Pslb,'--r')
   semilogy(10*log10(Es_No),Psub,'--r')
   axis([-20 20 1e-5 1])
   
   title(['Symbol Error Rate, M=' int2str(M)])
   xlabel('Es/N_o (dB)')
   ylabel('SER')
   hold off
   
subplot(1,2,2)
   semilogy(10*log10(Eb_No),BER,'oc')
   hold on
   semilogy(10*log10(Eb_No),Pbe,'-g')
   
   title(['Bit Error Rate, M=' int2str(M)])
   xlabel('Eb/N_o (dB)')
   ylabel('BER')
   hold off
   axis([-20 20 1e-5 1])

set(gcf,'Paperposition',[0.25 2.5 8 6])