numSymbs = 10000;
M = 4;

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

rolloff = 0.25;
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 = pskmod(data, M, 0, 'gray');
x = txFilter(modData, rolloff, span, sps);

%% Simulate chromatic dispersion
D = 20; % ps / (nm km)
lambda = 1550; % nm
z = 10; % km

xCD = chromaticDispersion(x, D, lambda, z, Tsamp);
xCD = normalizeEnergy(xCD, numSymbs, 1);


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

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

  yCDComp = CDCompensation(y, D, lambda, z, Tsamp);

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

  rSampled = r(sps*span/2+1:sps:(numSymbs + span/2) * sps);
  demodData = pskdemod(rSampled, M, 0, 'gray');

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

figure(1);
clf;

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

theoreticalPSK(EbN0_db, M, 'b', 'LineWidth', 1);
legend({'CD + AWGN + CD compensation', 'AWGN only'}, 'Location', 'southwest');

title(strcat(num2str(M), '-PSK with chromatic dispersion and compensation'));
grid on;
xlabel('$E_b/N_0$ (dB)');
ylabel('BER');

formatFigure;