No BSD License  

MPSK Simulation (BPSK,QPSK,8PSK)

by

 

20 Jun 2006 (Updated )

This m-file simulates MPSK (theoretical and simulated) with Gray coding

MPSK.m
%JC 6/20/06
%I was interested in the theory for 8PSK and came across this m-file written in 1999. It
%requires 2 function calls (graymapPSK and grayunmapPSK) which I have
%included and must be uncommented(single comments only) and loaded into your workspace.
%Remove the function calls from this m-file after you load into your
%workspace.
%It seems to work well and give valid results for BPSK and QPSK but I have a question
%about the setup for 8PSK and the solution for Pseint where Pse=Pseint.
%Also, the Gray coding seems different than what I'm used to.
%Some of you more theoretical minded may have some insight into the theory
%behind how Pseint was determined and comment. I don't have the
%communications toolbox-blockset and was wondering how it compared with Matlab's
%example simulation of 8PSK with Gray coding (in Simulink?)
%Each run takes about a minute and the loop iterations are shown in the
%command window.
%A  numerical example of a satellite link is shown using uncoded QPSK or 8PSK when the Rb is
%greater than the channel bandwidth Wc (Band-limited channel). It's a
%little wordy but I have tryed to provide enough information so you don't
%have to read between the lines.
% 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=symbols;
   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));%lower and upper bounds
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])
   
   %===============================================================
   %Satellite link numerical example 
   %===============================================================
   %Pr/No=EbRb/No=EsRs/NO
   %Wc=channel bandwidth
   %Rb=bit rate
   %Pr/No=power received/noise spectral density=dBW/dBW-Hz=dB-Hz
   %Remember that Pr/No is the level of the received signal plus noise at
   %the receiver's demodulator and includes the noise figure of the
   %receiver.
   %PB=bit error rate
   %PE=probability of error
   %Rs=symbol rate
   %Es=Energy of symbol
   %Eb=Energy of bit
   %log2(4)=2
   %log2(8)=3
   %===============================================================
   %(Wc=20MHz)  (Rb=30Mbits/sec)  (Pr/No=85dB-Hz)  (PB<1e-3 required)
   %Rb>Wc therefore band-limited channel requiring MPSK scheme(QPSK or 8PSK)
   %Therefore M=4 and Rs=Rb/(log2M)=30Mbits/sec/2=15Msymbols/sec which is <Wc=20MHz
   %Es/No=(log2M)Eb/No=log2M(Pr/NoRb)=2*(10exp8.5/30e6)=21.08=13.24dB
   %Eb/No=Pr/NoRb=(10exp8.5/30e6)=10.54=10*log10(10.54)=10.23dB
   %PE(M=4)=Q(sqrt(2*Es/No)*sin(pi/M))=Q(sqrt(2*21.08)*0.3826)=Q(2.484)
   %where Q(x)=.5*erfc(sqrt(x)/sqrt(2))
   %PE(M=4)=.5*erfc(1.5786/sqrt(2))=6.5e-3
   %PB=PE(4)/log2M=6.5e-3/2=3.25e-3 !!!!PB not low enough!!!!
   %Looks like QPSK didn't do the job so lets try 8PSK
   %Therefore M=8 and Rs=Rb/(log2M)=30Mbits/sec/3=10Msymbols/sec <Wc=20MHz
   %Es/No=(log2M)Eb/No=log2M(Pr/NoRb)=3*(10exp8.5/30e6)=31.62=15dB
   %Eb/No=Pr/NoRb=(10exp8.5/30e6)=10.54=10*log10(10.54)=10.23dB
   %PE(M=8)=2*Q(sqrt(2*Es/No)*sin(pi/M))=2*Q(sqrt(2*31.62)*0.3826)=2*Q(3.042)
   %PE(M=8)=2*.5*erfc(3.042/sqrt(2))=2.4e-3
   %PB=PE(8)/log2M=2.4e-3/3=8e-4 !!!!this PB is low enough and meets the required value!!!!
   %The power at the transmitter has been increased by 15dB-13.24dB=1.76dB
   %by using 8PSK over QPSK to meet the requirements of the system.
   
  %comment out only single comments
  % function symbols = graymapPSK(bits)

%K = size(bits,1);
%N = size(bits,2);

%switch K
%case 1             % BPSK
   %% maps 0=s1, 1=s0
   %% s0 = 1, s1 = -1
   %symbols = [bits*2-1;zeros(1,N)];
   %case 2
   %% maps 00=s0, 01=s1, 11=s2, 10=s3
   %% s0 = [1;0], s1 = [0;1], s2 = [-1;0], s3 = [0;-1]
%case 3
   %% maps 000=s0, 001=s1, 011=s2, 010=s3, 110=s4, 111=s5, 101=s6, 100=s7
   %% s0 = [1;0],  s1 = 1/sqrt(2)*[1;1],  s2 = [0;1],  s3 = 1/sqrt(2)*[-1;1],
   %% s4 = [-1;0], s5 = 1/sqrt(2)*[-1;-1],s6 = [0;-1], s7 = 1/sqrt(2)*[1;-1]
   %s = sum(bits,1);
   %s_even = [1-bits(1,:)-bits(2,:);bits(2,:)-bits(1,:)];
   %s_odd = (1/sqrt(2))*[-1+2*abs(bits(3,:)-bits(1,:));-1+2*abs(bits(3,:)-bits(2,:))];
   %symbols = s_even.*([1;1]*(s==0|s==2))+s_odd.*([1;1]*(s==1|s==3));
   %end
   
   %function bits = grayunmapPSK(symbols,M)

%K = log2(M);
%N = size(symbols,2);

%switch K
%case 1             % BPSK
   %% maps 0=s1, 1=s0
   %% s0 = 1, s1 = -1
   %bits = 0.5*(symbols(1,:)+1);
   %case 2
   %% maps 00=s0, 01=s1, 11=s2, 10=s3
   %% s0 = [1;0], s1 = [0;1], s2 = [-1;0], s3 = [0;-1]
   %bits = [0.5*(1-symbols(1,:)-symbols(2,:));0.5*(1+symbols(2,:)-symbols(1,:))];
   %case 3
   %% maps 000=s0, 001=s1, 011=s2, 010=s3, 110=s4, 111=s5, 101=s6, 100=s7
   %% s0 = [1;0],  s1 = 1/sqrt(2)*[1;1],  s2 = [0;1],  s3 = 1/sqrt(2)*[-1;1],
   %% s4 = [-1;0], s5 = 1/sqrt(2)*[-1;-1],s6 = [0;-1], s7 = 1/sqrt(2)*[1;-1]
   %bits_even = [0.5*(1-symbols(1,:)-symbols(2,:));0.5*(1+symbols(2,:)-symbols(1,:));...
         %abs(symbols(2,:))];
   %bits_odd = 0.5*[1-symbols(2,:)*sqrt(2);1-symbols(1,:)*sqrt(2);sqrt(2)*abs(symbols(1,:)+symbols(2,:))];
   %s = abs(symbols(1,:))>0.7 & abs(symbols(1,:))<0.8;
   %bits = bits_even.*([1;1;1]*(~s))+bits_odd.*([1;1;1]*(s));
   %end


   

   

Contact us