%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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          : 10th April 2010
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

clear
N = 10^6 % number of bits or symbols
M = 4; % QPSK constellation
k = log2(M); % number of bits per symbol

Eb_N0_dB = [0:15]; % multiple Eb/N0 values
Es_N0_dB = Eb_N0_dB + 10*log10(k);
 
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,2,N/2).';
    bintoDecConv = ones(N/2,1)*2.^[k-1:-1:0];	
    ipBitDec     = sum(ipBitReshape.*bintoDecConv,2);

    % converting to gray coded symbols
    ipBitGray    = bitxor(ipBitDec,floor(ipBitDec/2));
    ipPhaseGray    = 2*ipBitGray.'+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/4); 	

    % white gaussian noise, 0 mean
    n = 1/sqrt(2)*[randn(1,N/2) + j*randn(1,N/2)]; 
    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)*4/pi;
    quant_diffPhase = 2*floor(est_diffPhase/2)+1;  % quantizing
   
    % gray to binary
    quant_diffPhase(find(quant_diffPhase<0)) = quant_diffPhase(find(quant_diffPhase<0)) + 8;
    bin_diffPhase     = floor(bitxor(quant_diffPhase,floor(quant_diffPhase/2))/2);
    estBit_noncoh     = (dec2bin(bin_diffPhase.')).'; 
    estBit_noncoh     = str2num(estBit_noncoh(1:end).').';

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

    % theoretical BER computation 
    a = sqrt(2*10.^(Eb_N0_dB(ii)/10)*(1-sqrt(1/2)));
    b = sqrt(2*10.^(Eb_N0_dB(ii)/10)*(1+sqrt(1/2)));
    
    k_bessel = 0:10; 
    temp = exp(-((a.^2+b.^2)/2)).*sum((a/b).^k_bessel.*besseli(k_bessel,a*b));
    theoryBer_dqpsk_noncoh(ii) = temp - 0.5*besseli(0,a*b)*exp(-((a.^2+b.^2)/2));

end

simBer_dqpsk_noncoh = nErr_dqpsk_noncoh/N;

close all
figure
semilogy(Eb_N0_dB,theoryBer_dqpsk_noncoh,'bx-');
hold on
semilogy(Eb_N0_dB,simBer_dqpsk_noncoh,'ms-');
axis([0 14 10^-5 0.5])
grid on
legend('theory-dqpsk-noncoh', 'simulation-dqpsk-noncoh');
xlabel('Eb/No, dB')
ylabel('Bit Error Rate')
title('Bit error probability curve for non coherent demodulation of pi/4 DQPSK')