4 Rsym = 2.5e10; % symbol rate (sym/sec)
5 Tsym = 1 / Rsym; % symbol period (sec)
8 span = 6; % filter span
9 sps = 2; % samples per symbol
11 fs = Rsym * sps; % sampling freq (Hz)
14 t = (0 : 1 / fs : numSymbs / Rsym + (1.5 * span * sps - 1) / fs).';
17 power_dBm = -3:0.25:4;
18 power = 10 .^ (power_dBm / 10) * 1e-3; % watts
20 Es = power * Tsym; % joules
21 Eb = Es / log2(M); % joules
23 N0ref_dB = 10; % Eb/N0 at power = 1mW
24 %% Fix N0, such that Eb/N0 = N0ref_dB at power = 1mW
25 N0 = 1e-3 * Tsym / (log2(M) * 10 ^ (N0ref_dB / 10)); % joules
28 plotlen = length(power);
30 ber = zeros(1, plotlen);
32 data = randi([0 M - 1], numSymbs, 1);
33 modData = pskmod(data, M, pi / M, 'gray');
36 %% Chromatic dispersion
37 D = 17; % ps / (nm km)
43 snr = Es(i) / sps / N0;
44 snr_dB = 10 * log10(snr);
46 x = txFilter(modData, rolloff, span, sps);
47 %% Now, sum(abs(x) .^ 2) / length(x) should be 1.
48 %% We can set its power simply by multiplying.
49 x = sqrt(power(i)) * x;
51 %% We can now do split-step Fourier.
52 gamma = 1.2; % watt^-1 / km
53 %%stepnum = round(40 * z * gamma); % Nonlinear Fiber optics, App B
55 xCD = splitstepfourier(x, D, lambda, z, Tsamp, gamma, stepnum);
57 y = awgn(xCD, snr, power(i), 'linear');
60 r = rxFilter(y, rolloff, span, sps);
61 rCDComp = CDCompensation(r, D, lambda, z, Tsamp);
62 rCDComp = normalizeEnergy(rCDComp, numSymbs*sps, 1);
64 rSampled = rCDComp(sps*span/2+1:sps:(numSymbs+span/2)*sps);
65 rNoCompSampled = r(sps*span/2+1:sps:(numSymbs+span/2)*sps);
67 %% rotate rNoCompSampled to match original data
68 theta = angle(-sum(rNoCompSampled .^ M)) / M;
69 %% if theta approx +pi/M, wrap to -pi/M
70 if abs(theta - pi / M) / (pi / M) < 0.1
73 rNoCompSampled = rNoCompSampled .* exp(-j * theta);
76 %% Not entirely sure why, but after using FFT instead of time-domain
77 %% convolution for simulating CD, we now need to do the same rotation
78 %% for rSampled as well, but this time with a positive rotation.
79 theta = angle(-sum(rSampled .^ M)) / M;
80 if abs(theta + pi / M) / (pi / M) < 0.1
83 rSampled = rSampled .* exp(-1j * theta);
87 adaptFilterOut = adaptiveCMA(rSampled);
89 demodAdapt = pskdemod(adaptFilterOut, M, pi / M, 'gray');
90 [~, ber(i)] = biterr(data, demodAdapt);
96 %% Plot simulated results
97 semilogy(power_dBm, ber, 'Color', [0, 0.6, 0], 'LineWidth', 2);
100 title({'CD + Kerr + CD compensation', ...
101 strcat(['$D = 17$ ps/(nm km), $z = ', num2str(z), '$ km']), ...
102 strcat(['$E_b/N_0 = ', num2str(N0ref_dB), '$ dB at 1 mW'])});
104 xlabel('Optical power (dBm)');