%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% All rights reserved by Krishna Pillai, http://www.dsplog.com
% The file may not be re-distributed without explicit authorization
% from Krishna Pillai.
% Checked for proper operation with Octave Version 3.0.0
% Author        : Krishna Pillai
% Email         : krishna@dsplog.com
% Version       : 1.0
% Date          : 19 July 2010
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Simple Matlab/Octave code for non-coherent demodulation of
% differentially encoded 8 phase shift keying (D8PSK)

clear
N = 9*10^5 % number of bits or symbols
M = 8; % QPSK constellation
k = log2(M); % number of bits per symbol

Eb_N0_dB = [0:20]; % multiple Eb/N0 values
Es_N0_dB = Eb_N0_dB + 10*log10(k);
phaseArray   = [1 3  15 13 7 5 9 11];
bitDecArray  = [0 1 5 4 6 7 3 2 ];
 
for ii = 1:length(Eb_N0_dB)

    % generating random binary sequence	
    ip = rand(1,N)>0.5; % generating 0,1 with equal probability

    % grouping to form QPSK symbols	
    ipBitReshape = reshape(ip,3,N/3).';
    bintoDecConv = ones(N/3,1)*2.^[k-1:-1:0];	
    ipBitDec     = sum(ipBitReshape.*bintoDecConv,2);

    % converting to gray coded symbols
    ipPhaseGray    = phaseArray(ipBitDec'+1);

    % generating differential modulated symbols
    % phi[k] = phi[k-1] + Dphi[k]
    diffPhase = filter([ 1 ],[1 -1],ipPhaseGray); % start with 0 phase 
    dqpsk_symbols = exp(j*diffPhase*pi/8); 	

    % white gaussian noise, 0 mean
    n = 1/sqrt(2)*[randn(1,N/3) + j*randn(1,N/3)]; 
    y_dqpsk = dqpsk_symbols + 10^(-(Es_N0_dB(ii))/20)*n; % additive white gaussian noise

    % non-coherent demodulation
    estPhase  = angle(y_dqpsk);
    % Dphi[k] = phi[k] - phi[k-1]
    est_diffPhase = filter([1 -1],1,estPhase)*8/pi;
    quant_diffPhase = 2*floor(est_diffPhase/2)+1;  % quantizing
    quant_diffPhase(find(quant_diffPhase<0)) = quant_diffPhase(find(quant_diffPhase<0)) + 16;
   
    % gray to binary
    estBitDec = bitDecArray(floor(quant_diffPhase/2)+1); 		
    estBit_noncoh     = (dec2bin(estBitDec.')).'; 
    estBit_noncoh     = str2num(estBit_noncoh(1:end).').';

    % counting errors 
    nErr_dqpsk_noncoh(ii) = size(find([ip - estBit_noncoh]),2); % counting the number of errors

end

simBer_dqpsk_noncoh = nErr_dqpsk_noncoh/N;

close all
figure
semilogy(Eb_N0_dB,simBer_dqpsk_noncoh,'ms-');
axis([0 20 10^-5 0.5])
grid on
legend('simulation-d8psk-noncoh');
xlabel('Eb/No, dB')
ylabel('Bit Error Rate')
title('Bit error probability curve for non coherent demodulation of pi/8 D8PSK')