%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 0.7 & abs(symbols(1,:))<0.8;
%bits = bits_even.*([1;1;1]*(~s))+bits_odd.*([1;1;1]*(s));
%end