%          ----------------------------------------
%  +++++    Simulation of DPSK Demodulator     +++++++
%  For simplicity only the detector (after sampling) is simulated.
%          ----------------------------------------   
clear all ; close all;

E = 1 ;       % normalized to unit energy
snr_dB = 18 ;

N = 10000 ;    % no of bits  ...
M = 8 ;        % no. of levels ...
phi = 12*pi/180  ;     % off-set in radians ..

snr = exp(snr_dB*log(10)/10) ;
sig = E/sqrt(2*snr) ;       % this is detector noise std deviation

% the detector signals and symbol error ....
err = 0 ;
const = 0.0 ;
r0(1) = cos(const) + sig*randn ;  % initial values 
r1(1) = sin(const) + sig*randn ;
for i=2:N+1
   tt = rand;            %  tt is uniformly distributed in 0->1 .
   src = floor(rand*M) ; %  simulate symbols 0 -- M-1 ....
   const = const + 2*pi*src/M ;  % find the change in angle ...
   if (const > 2*pi)
       const = const - 2*pi ;
   end ;
   r0(i) = cos(const+phi) + sig*randn ;  % WGN 
   r1(i) = sin(const+phi) + sig*randn ;
   difr(i) = (r0(i) + j*r1(i))*(r0(i-1) - j*r1(i-1)) ;
   ang = atan2(r1(i),r0(i)) - atan2(r1(i-1),r0(i-1)) ;
   ang = atan2(imag(difr(i)),real(difr(i))) ;
   if(ang < 0)
       ang = ang + 2*pi ;   % angle in 0 --> 2*pi
   end ;
   det = round(M*ang/(2*pi)) ;
   if(det == M ) 
       det = 0 ;
   end ;
   if(det ~= src) 
      err = err + 1 ;
   end ;
end;
SER = err/N ;

% calculte therotical BER ...in a SNR (db) range
dbR = 0:1:20 ;
for i = 1:length(dbR)
   snr = exp(dbR(i)*log(10)/10) ;
   SERTpsk(i) = erfc(sqrt(snr)*sin(pi/M)) ;
   SERTdpsk(i) = erfc(sqrt(snr/2)*sin(pi/M)) ;
end ;

figure(1); plot(r0,r1,'.');
hold on ; 
for k=1:M
 plot([0 2*cos((2*k-1)*pi/M)], [0 2*sin((2*k-1)*pi/M)],'r--') ; 
end ;
hold off ;
axis([-1.5 1.5 -1.5 1.5]); 
grid;
title('Received signal constellations');
xlabel('real part');
ylabel('imaginary part');

figure(2); plot(real(difr),imag(difr),'.');
hold on ; 
for k=1:M
 plot([0 2*cos((2*k-1)*pi/M)], [0 2*sin((2*k-1)*pi/M)],'r--') ; 
end ;
hold off ;
axis([-1.5 1.5 -1.5 1.5]); 
grid;
title('Differential signal constellations');
xlabel('real part');
ylabel('imaginary part');

figure(3); semilogy(dbR, SERTdpsk, dbR, SERTpsk, '--', snr_dB, SER,'r*') ; 
grid ;
xlabel('SNR in dB');
ylabel('Symbol Error Probability, Pm');
title('Theoretical and simulated SER for DPSK');
legend('DPSK Thero','PSK Thero','DPSK Simo');