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

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

span = 6; % filter span
sps = 4; % 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);

ber = zeros(1, plotlen);

data = randi([0 M - 1], numSymbs, 1);
%%modData = dpskmod(data, M, pi / M, 'gray');
modData = dpskmod(data, M, 0, 'gray');

x = txFilter(modData, rolloff, span, sps);

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


avgSa = 40;


[xPN, pTxLO] = phaseNoise(x, linewidthTx, linewidthLO, Tsamp);

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

  y = awgn(xPN, snr, 'measured');

  r = rxFilter(y, rolloff, span, sps);
  %% normalize energy
  r = normalizeEnergy(r, numSymbs, 1 + noiseEnergy);

  rSampled = r(sps*span/2+1:sps:(numSymbs+span/2)*sps);
  [rSaPhEq, phiests] = phaseNoiseCorr(rSampled, M, avgSa);

  adaptiveFilterOut = adaptiveCMA(rSaPhEq.');

  demodData = dpskdemod(rSaPhEq, M, 0, 'gray').';
  demodAdapt = dpskdemod(adaptiveFilterOut, M, 0, 'gray');

  [~, ber(i)] = biterr(data, demodData);

  if EbN0_db(i) == 8
    figure(1234);
    plot(repelem(-phiests, sps));
    hold on;
    plot(pTxLO);
    legend('estimate', 'actual');
    hold off;

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


figure(1);
clf;

%% Plot simulated results
semilogy(EbN0_db, ber, 'r', 'LineWidth', 2);
hold on;

theoreticalPSK(EbN0_db, M, 'b', 'LineWidth', 1);
legend({'Simulated phase noise + correction', 'Theoretical AWGN'}, ...
       'Location', 'southwest');


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

formatFigure;