| 1 | numSymbs = 2^16; |
| 2 | M = 4; |
| 3 | |
| 4 | Rsym = 2.5e10; % symbol rate (sym/sec) |
| 5 | Tsym = 1 / Rsym; % symbol period (sec) |
| 6 | |
| 7 | rolloff = 0.25; |
| 8 | span = 6; % filter span |
| 9 | sps = 8; % samples per symbol |
| 10 | |
| 11 | fs = Rsym * sps; % sampling freq (Hz) |
| 12 | Tsamp = 1 / fs; |
| 13 | |
| 14 | t = (0 : 1 / fs : numSymbs / Rsym + (1.5 * span * sps - 1) / fs).'; |
| 15 | |
| 16 | power_dBm = -6:1:4; |
| 17 | %%power_dBm = 0; |
| 18 | power = 10 .^ (power_dBm / 10) * 1e-3; % watts |
| 19 | |
| 20 | Es = power * Tsym; % joules |
| 21 | Eb = Es / log2(M); % joules |
| 22 | |
| 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 |
| 27 | |
| 28 | plotlen = length(power); |
| 29 | |
| 30 | ber = zeros(1, plotlen); |
| 31 | |
| 32 | data = randi([0 M - 1], numSymbs, 1); |
| 33 | %%modData = dpskmod(data, M, 0, 'gray'); |
| 34 | modData = pskmod(data, M, 0, 'gray'); |
| 35 | for i = 2:numSymbs |
| 36 | modData(i) = modData(i) * modData(i-1); |
| 37 | end |
| 38 | |
| 39 | |
| 40 | %% Chromatic dispersion |
| 41 | D = 17; % ps / (nm km) |
| 42 | lambda = 1550; % nm |
| 43 | z = 100; % km |
| 44 | |
| 45 | |
| 46 | linewidthTx = 0; % Hz |
| 47 | linewidthLO = 1e6; % Hz |
| 48 | |
| 49 | |
| 50 | TsampOrig = Tsamp; |
| 51 | |
| 52 | x_P1 = txFilter(modData, rolloff, span, sps); |
| 53 | |
| 54 | |
| 55 | for i = 1:plotlen |
| 56 | sps = 8; |
| 57 | Tsamp = TsampOrig; |
| 58 | |
| 59 | snr = Es(i) / sps / N0; |
| 60 | snr_dB = 10 * log10(snr); |
| 61 | |
| 62 | %%x = txFilter(modData, rolloff, span, sps); |
| 63 | %% Now, sum(abs(x) .^ 2) / length(x) should be 1. |
| 64 | %% We can set its power simply by multiplying. |
| 65 | x = sqrt(power(i)) * x_P1; |
| 66 | |
| 67 | %% We can now do split-step Fourier. |
| 68 | gamma = 1.2; % watt^-1 / km |
| 69 | |
| 70 | |
| 71 | xCDKerr = splitstepfourier(x, D, lambda, z, Tsamp, gamma); |
| 72 | |
| 73 | xpn = phaseNoise(xCDKerr, linewidthTx, linewidthLO, Tsamp); |
| 74 | |
| 75 | y = awgn(xpn, snr_dB, 'measured', 'db'); |
| 76 | %y = xCDKerr; |
| 77 | |
| 78 | r = rxFilter(y, rolloff, span, sps); |
| 79 | sps = 2; |
| 80 | Tsamp = Tsamp * 4; |
| 81 | |
| 82 | rCDComp = CDCompensation(r, D, lambda, z, Tsamp); |
| 83 | rCDComp = normalizeEnergy(rCDComp, numSymbs * sps, 1); |
| 84 | |
| 85 | rSampled = rCDComp(2:2:end); |
| 86 | |
| 87 | %% adaptive filter |
| 88 | [adaptFilterOut, convergeIdx] = adaptiveCMA(rSampled); |
| 89 | |
| 90 | pncorr = phaseNoiseCorr(adaptFilterOut, M, 0, 40).'; |
| 91 | |
| 92 | demodAdapt = pskdemod(pncorr, M, 0, 'gray'); |
| 93 | remod = pskmod(demodAdapt, M, 0, 'gray'); |
| 94 | delayed = [1; remod(1:end-1)]; |
| 95 | demod = pskdemod(remod .* conj(delayed), M, 0, 'gray'); |
| 96 | |
| 97 | if convergeIdx < Inf |
| 98 | [~, ber(i)] = biterr(data(convergeIdx:end), demod(convergeIdx:end)); |
| 99 | else |
| 100 | [~, ber(i)] = biterr... |
| 101 | (data(ceil(0.8*numSymbs):end), ... |
| 102 | demod(ceil(0.8*numSymbs):end)); |
| 103 | end |
| 104 | end |
| 105 | |
| 106 | ber |
| 107 | |
| 108 | |
| 109 | figure(1); |
| 110 | clf; |
| 111 | |
| 112 | %% Plot simulated results |
| 113 | qp = 20 * log10(erfcinv(2*ber)*sqrt(2)); |
| 114 | plot(power_dBm, qp, 'Color', [0, 0.6, 0], 'LineWidth', 2); |
| 115 | hold on; |
| 116 | |
| 117 | title({'CD + Kerr + CD compensation', ... |
| 118 | strcat(['$D = 17$ ps/(nm km), $z = ', num2str(z), '$ km'])}); |
| 119 | grid on; |
| 120 | xlabel('Optical power (dBm)'); |
| 121 | ylabel('$20 \log_{10}\left(\sqrt{2}\mathrm{erfc}^{-1}(2 BER)\right)$'); |
| 122 | |
| 123 | formatFigure; |