numSymbs = 5e4;
M = 4;
rolloff = 0.5;

Rsym = 2.5e10; % symbol rate (sym/sec)

span = 6; % filter span
sps = 8; % samples per symbol

fs = Rsym * sps; % sampling freq (Hz)
Tsamp = 1 / fs;

t = (0 : 1 / fs : numSymbs / Rsym + (1.5 * span * sps - 1) / fs).';


EbN0_db = 0:0.2:14;
EbN0 = 10 .^ (EbN0_db ./ 10);

Es = 1;
Eb = Es / log2(M);
N0 = Eb ./ EbN0;

EsN0 = EbN0 .* log2(M);
EsN0_db = 10 .* log10(EsN0);

plotlen = length(EbN0);

berPSK = zeros(1, plotlen);
berDEPSK = zeros(1, plotlen);
berDPSK = zeros(1, plotlen);

data = randi([0 M - 1], numSymbs, 1);

pskSym = pskmod(data, M, pi / M, 'gray');
%% DEPSK: Part VII, M.G. Taylor (2009)
depskSym = pskmod(data, M, 0, 'gray');
for i = 2:numSymbs
  depskSym(i) = depskSym(i) * depskSym(i-1);
end

dpskSym = dpskmod(data, M, pi / M, 'gray');

xPSK = txFilter(pskSym, rolloff, span, sps);
xDEPSK = txFilter(depskSym, rolloff, span, sps);
xDPSK = txFilter(dpskSym, rolloff, span, sps);

linewidthTx = 0; % Hz
linewidthLO = 5e6; % Hz
%linewidthLO = Rsym * 1e-3;

iterations = 1;
avgSa = 40;

TsampOrig = Tsamp;

for it = 1 : iterations
  [xPSKpn, pTxLoPSK] = phaseNoise(xPSK, linewidthTx, linewidthLO, Tsamp);
  [xDEPSKpn, pTxLoDEPSK] = phaseNoise(xDEPSK, linewidthTx, linewidthLO, Tsamp);
  [xDPSKpn, pTxLoDPSK] = phaseNoise(xDPSK, linewidthTx, linewidthLO, Tsamp);

  for i = 1:plotlen
    Tsamp = TsampOrig;
    sps = 8;

    snr = EbN0_db(i) + 10 * log10(log2(M)) - 10 * log10(sps);
    noiseEnergy = 10 ^ (-snr / 10);

    yPSK = awgn(xPSKpn, snr, 'measured');
    yDEPSK = awgn(xDEPSKpn, snr, 'measured');
    yDPSK = awgn(xDPSKpn, snr, 'measured');

    rPSK = rxFilter(yPSK, rolloff, span, sps);
    rDEPSK = rxFilter(yDEPSK, rolloff, span, sps);
    rDPSK = rxFilter(yDPSK, rolloff, span, sps);

    sps = 2;
    Tsamp = TsampOrig * 4;

    rPSKSamp = rPSK(1:2:end);
    rDEPSKSamp = rDEPSK(1:2:end);
    rDPSKSamp = rDPSK(1:2:end);

    [rPSKSampEq, phiestsPSK] = phaseNoiseCorr(rPSKSamp, M, pi/M, avgSa);
    [rDEPSKSampEq, phiestsDEPSK] = phaseNoiseCorr(rDEPSKSamp, M, 0, avgSa);

    demodPSK = pskdemod(rPSKSampEq, M, pi/M, 'gray').';
    %% The decoding method described in Taylor (2009)
    %% works on the complex symbols, i.e. after taking
    %% the nearest symbol in the constellation, but before
    %% converting them back to integers/bits.
    %% MATLAB's pskdemod() does not provide this intermediate
    %% result, so to be lazy, a pskmod() call is performed
    %% to obtain the complex symbols.
    demodDEPSK = pskdemod(rDEPSKSampEq, M, 0, 'gray').';
    remodDEPSK = pskmod(demodDEPSK, M, 0, 'gray');
    delayed = [1; remodDEPSK(1:end-1)];
    demodDEPSK = pskdemod(remodDEPSK .* conj(delayed), M, 0, 'gray');

    demodDPSK = dpskdemod(rDPSKSamp, M, pi/M, 'gray');

    [~, ber] = biterr(data, demodPSK);
    berPSK(i) = berPSK(i) + ber / iterations;
    [~, ber] = biterr(data, demodDEPSK);
    berDEPSK(i) = berDEPSK(i) + ber / iterations;
    [~, ber] = biterr(data, demodDPSK);
    berDPSK(i) = berDPSK(i) + ber / iterations;

    if EbN0_db(i) == 8 && it == 1
      figure(1234);
      plot(repelem(-phiestsPSK, 8));
      hold on;
      plot(pTxLoPSK);
      legend('estimate', 'actual');
      hold off;

      figure(1);
      scatterplot(rPSKSampEq);
      title('rPSKSampEq');
    end

  end
end


figure(1);
clf;

%% Plot simulated results
semilogy(EbN0_db, berPSK, 'r', 'LineWidth', 1.5);
hold on;
semilogy(EbN0_db, berDEPSK, 'c', 'LineWidth', 2);
semilogy(EbN0_db, berDPSK, 'Color', [0, 0.6, 0], 'LineWidth', 2.5);

theoreticalPSK(EbN0_db, M, 'b', 'LineWidth', 1);
DEPSKTheoretical = berawgn(EbN0_db, 'psk', M, 'diff');
semilogy(EbN0_db, DEPSKTheoretical, 'Color', [1, 0.6, 0], 'LineWidth', 1);
DPSKTheoretical = berawgn(EbN0_db, 'dpsk', M);
semilogy(EbN0_db, DPSKTheoretical, 'm', 'LineWidth', 1);

legend({'PSK with Viterbi-Viterbi', ...
        'DEPSK with Viterbi-Viterbi', ...
        'DPSK', ...
        'Theoretical PSK over AWGN', ...
        'Theoretical DEPSK over AWGN', ...
        'Theoretical DPSK over AWGN'}, ...
       'Location', 'southwest');

title({'QPSK with phase nosie and correction', ...
       strcat('$10^{', num2str(log10(numSymbs * log2(M))), ...
              '}$~bits, LO~', ...
              num2str(linewidthLO / 1e6), '~MHz, blocksize~', ...
              num2str(avgSa), '~Sa')});
grid on;
xlabel('$E_b/N_0$ (dB)');
ylabel('BER');

formatFigure;