Commit | Line | Data |
---|---|---|
5fae0077 AIL |
1 | numSymbs = 2^14; |
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/2 / Rsym - 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 | linewidthTx = 0; % Hz | |
46 | linewidthLO = 1e6; % Hz | |
47 | ||
48 | TsampOrig = Tsamp; | |
49 | ||
50 | sig_x = txFilter(modData(1:numSymbs/2), rolloff, span, sps); | |
51 | sig_y = txFilter(modData(numSymbs/2+1:end), rolloff, span, sps); | |
52 | ||
53 | for i = 1:plotlen | |
54 | sps = 8; | |
55 | Tsamp = TsampOrig; | |
56 | ||
57 | snr = Es(i) / sps / N0; | |
58 | snr_dB = 10 * log10(snr); | |
59 | ||
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)); | |
66 | ||
67 | rot_omega = 1e3; % rad/s | |
68 | rot_phi = 2; | |
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); | |
73 | ||
74 | %% We can now do split-step Fourier. | |
75 | gamma = 1.2; % watt^-1 / km | |
76 | ||
77 | [xCDKerr, yCDKerr] = ssf_pdm(rot_x, rot_y, ... | |
78 | D, lambda, z, Tsamp, gamma); | |
79 | ||
80 | xpn = phaseNoise(xCDKerr, linewidthTx, linewidthLO, Tsamp); | |
81 | ypn = phaseNoise(yCDKerr, linewidthTx, linewidthLO, Tsamp); | |
82 | ||
83 | xout = awgn(xpn, snr_dB, 'measured', 'db'); | |
84 | yout = awgn(ypn, snr_dB, 'measured', 'db'); | |
85 | ||
86 | rx = rxFilter(xout, rolloff, span, sps); | |
87 | ry = rxFilter(yout, rolloff, span, sps); | |
88 | sps = 2; | |
89 | Tsamp = Tsamp * 4; | |
90 | ||
91 | rxCDComp = CDCompensation(rx, D, lambda, z, Tsamp); | |
92 | ryCDComp = CDCompensation(ry, D, lambda, z, Tsamp); | |
93 | ||
94 | rxSampled = rxCDComp(1:2:end); | |
95 | rySampled = ryCDComp(1:2:end); | |
96 | ||
97 | %% adaptive filter | |
98 | [xCMA, yCMA] = pdm_adaptiveCMA(rxSampled, rySampled); | |
99 | ||
100 | xpncorr = phaseNoiseCorr(xCMA, M, 0, 40).'; | |
101 | ypncorr = phaseNoiseCorr(yCMA, M, 0, 40).'; | |
102 | ||
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'); | |
107 | clear remodx | |
108 | clear delayx | |
109 | ||
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'); | |
114 | clear remody | |
115 | clear delayy | |
116 | ||
117 | ||
118 | [~, ber(i)] = biterr(data, [demodx; demody]); | |
119 | end | |
120 | ||
121 | ber | |
122 | ||
123 | figure; | |
124 | clf; | |
125 | ||
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); | |
129 | hold on; | |
130 | ||
131 | title({'CD + Kerr + CD compensation', ... | |
132 | strcat(['$D = 17$ ps/(nm km), $z = ', num2str(z), '$ km'])}); | |
133 | grid on; | |
134 | xlabel('Optical power (dBm)'); | |
135 | ylabel('$20 \log_{10}\left(\sqrt{2}\mathrm{erfc}^{-1}(2 BER)\right)$'); | |
136 | ||
137 | formatFigure; |