numSymbs = 10000;
M = 4;

Rsym = 2.5e10; % symbol rate (sym/sec)
rolloff = 0.5;
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 = 8;
EbN0 = 10 .^ (EbN0_db ./ 10);

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

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


data = randi([0 M - 1], numSymbs, 1);
pskSym = pskmod(data, M, pi/M, 'gray');
dpskSym = dpskmod(data, M, pi/M, 'gray');

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

linewidthTx = 0; % Hz
linewidthLO = 10e6; % Hz

[xPSKpn, pTxLoPSK] = phaseNoise(xPSK, linewidthTx, linewidthLO, Tsamp);
[xDPSKpn, pTxLoDPSK] = phaseNoise(xDPSK, linewidthTx, linewidthLO, Tsamp);


snr = EbN0_db + 10 * log10(log2(M)) - 10 * log10(sps);

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

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

sps = 2;
Tsamp = Tsamp * 4;

rPSKSa = rPSK(1:2:end);
rDPSKSa = rDPSK(1:2:end);

[rPSKSaPhEq, phiestsPSK] = phaseNoiseCorr(rPSKSa, M, pi/M, 80);

demodPSK = pskdemod(rPSKSaPhEq, M, pi/M, 'gray');
demodDPSK = dpskdemod(rDPSKSa, M, pi/M, 'gray');


[bitErrors, ber] = biterr(data, demodPSK.')


figure(2);
plot(t(1:80000), repelem(-phiestsPSK, 8));
hold on;
plot(t(1:80000), pTxLoPSK(1:80000));
legend('estimate', 'actual');
title('Phase noise estimation');
hold off;
return
figure(3);
plot(t(1:length(x)), real(normalizeEnergy(x, numSymbs*sps, 1)));
hold on;
%%plot(t, real(rPhaseEq), 'r');
sampledTimes = t(sps*span/2+1:sps:(numSymbs+span/2)*sps);
%%plot(sampledTimes, real(rSampled), 'x');
plot(sampledTimes, real(normalizeEnergy(rSaPhEq, numSymbs, 1)), 'x');
legend('original signal', 'corrected received signal');
title('Phase noise correction, linewidth 1 MHz, $E_b/N_0=8$ dB');
ylabel('Real part of signals');
axis([t(1) t(300) -Inf +Inf]);
hold off;
formatFigure;