4 Rsym = 2.5e10; % symbol rate (sym/sec)
5 Tsym = 1 / Rsym; % symbol period (sec)
8 span = 6; % filter span
9 sps = 8; % samples per symbol
11 fs = Rsym * sps; % sampling freq (Hz)
14 t = (0 : 1 / fs : numSymbs/2 / Rsym - 1/fs).';
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
26 %% At current settings, N0 = 0.002 pJ
28 plotlen = length(power);
30 ber = zeros(1, plotlen);
32 data = randi([0 M - 1], numSymbs, 1);
33 %%modData = dpskmod(data, M, 0, 'gray');
34 modData = pskmod(data, M, 0, 'gray');
36 modData(i) = modData(i) * modData(i-1);
40 %% Chromatic dispersion
41 D = 17; % ps / (nm km)
46 linewidthLO = 1e6; % Hz
50 sig_x = txFilter(modData(1:numSymbs/2), rolloff, span, sps);
51 sig_y = txFilter(modData(numSymbs/2+1:end), rolloff, span, sps);
57 snr = Es(i) / sps / N0;
58 snr_dB = 10 * log10(snr);
60 %%x = txFilter(modData, rolloff, span, sps);
61 %% Now, sum(abs(x) .^ 2) / length(x) should be 1.
62 %% We can set its power simply by multiplying.
63 %%x = sqrt(power(i)) * x;
64 txx = sig_x * sqrt(power(i));
65 txy = sig_y * sqrt(power(i));
67 rot_omega = 1e3; % rad/s
69 rot_x = txx .* cos(rot_omega * t) + ...
70 txy .* sin(rot_omega * t) * exp(-1j * rot_phi);
71 rot_y = txx .* -sin(rot_omega * t) * exp(1j * rot_phi) + ...
72 txy .* cos(rot_omega * t);
74 %% We can now do split-step Fourier.
75 gamma = 1.2; % watt^-1 / km
77 [xCDKerr, yCDKerr] = ssf_pdm(rot_x, rot_y, ...
78 D, lambda, z, Tsamp, gamma);
80 xpn = phaseNoise(xCDKerr, linewidthTx, linewidthLO, Tsamp);
81 ypn = phaseNoise(yCDKerr, linewidthTx, linewidthLO, Tsamp);
83 xout = awgn(xpn, snr_dB, 'measured', 'db');
84 yout = awgn(ypn, snr_dB, 'measured', 'db');
86 rx = rxFilter(xout, rolloff, span, sps);
87 ry = rxFilter(yout, rolloff, span, sps);
91 rxCDComp = CDCompensation(rx, D, lambda, z, Tsamp);
92 ryCDComp = CDCompensation(ry, D, lambda, z, Tsamp);
94 rxSampled = rxCDComp(1:2:end);
95 rySampled = ryCDComp(1:2:end);
98 [xCMA, yCMA] = pdm_adaptiveCMA(rxSampled, rySampled);
100 xpncorr = phaseNoiseCorr(xCMA, M, 0, 40).';
101 ypncorr = phaseNoiseCorr(yCMA, M, 0, 40).';
103 demodx = pskdemod(xpncorr, M, 0, 'gray');
104 remodx = pskmod(demodx, M, 0, 'gray');
105 delayx = [1; remodx(1:end-1)];
106 demodx = pskdemod(remodx .* conj(delayx), M, 0, 'gray');
110 demody = pskdemod(ypncorr, M, 0, 'gray');
111 remody = pskmod(demody, M, 0, 'gray');
112 delayy = [1; remody(1:end-1)];
113 demody = pskdemod(remody .* conj(delayy), M, 0, 'gray');
118 [~, ber(i)] = biterr(data, [demodx; demody]);
126 %% Plot simulated results
127 qp = 20 * log10(erfcinv(2*ber)*sqrt(2));
128 plot(power_dBm, qp, 'Color', [0, 0.6, 0], 'LineWidth', 2);
131 title({'CD + Kerr + CD compensation', ...
132 strcat(['$D = 17$ ps/(nm km), $z = ', num2str(z), '$ km'])});
134 xlabel('Optical power (dBm)');
135 ylabel('$20 \log_{10}\left(\sqrt{2}\mathrm{erfc}^{-1}(2 BER)\right)$');