Commit | Line | Data |
---|---|---|
f9a73e9e AIL |
1 | numSymbs = 5e5; |
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 = 2; % 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 | ||
17 | %%power_dBm = -3:0.2:4; | |
18 | power_dBm = [0]; | |
19 | power = 10 .^ (power_dBm / 10) * 1e-3; % watts | |
20 | ||
21 | Es = power * Tsym; % joules | |
22 | Eb = Es / log2(M); % joules | |
23 | ||
24 | N0ref_db = 10; % Eb/N0 at power = 1mW | |
25 | %% Fix N0, such that Eb/N0 = N0ref_db at power = 1mW | |
26 | N0 = 1e-3 * Tsym / (log2(M) * 10 ^ (N0ref_db / 10)); % joules | |
27 | ||
28 | ||
29 | plotlen = length(power); | |
30 | ||
31 | ber = zeros(1, plotlen); | |
32 | ||
33 | data = randi([0 M - 1], numSymbs, 1); | |
34 | modData = pskmod(data, M, pi / M, 'gray'); | |
35 | ||
36 | ||
37 | %% Chromatic dispersion | |
38 | D = 17; % ps / (nm km) | |
39 | lambda = 1550; % nm | |
40 | z = 600; % km | |
41 | ||
42 | ||
43 | for i = 1:plotlen | |
44 | snr = Es(i) / sps / N0; | |
45 | snr_dB = 10 * log10(snr); | |
46 | ||
47 | x = txFilter(modData, rolloff, span, sps); | |
48 | %% Now, sum(abs(x) .^ 2) / length(x) should be 1. | |
49 | %% We can set its power simply by multiplying. | |
50 | x = sqrt(power(i)) * x; | |
51 | ||
52 | %% We can now do split-step Fourier. | |
53 | gamma = 1.2; % watt^-1 / km | |
54 | %%stepnum = round(40 * z * gamma); % Nonlinear Fiber optics, App B | |
55 | stepnum = 100; | |
56 | xCD = splitstepfourier(x, D, lambda, z, Tsamp, gamma, stepnum); | |
57 | ||
58 | y = awgn(xCD, snr, power(i), 'linear'); | |
59 | %%y = xCD; | |
60 | ||
61 | r = rxFilter(y, rolloff, span, sps); | |
62 | rCDComp = CDCompensation(r, D, lambda, z, Tsamp); | |
63 | rCDComp = normalizeEnergy(rCDComp, numSymbs*sps, 1); | |
64 | ||
65 | rSampled = rCDComp(sps*span/2+1:sps:(numSymbs+span/2)*sps); | |
66 | rNoCompSampled = r(sps*span/2+1:sps:(numSymbs+span/2)*sps); | |
67 | ||
68 | %% rotate rNoCompSampled to match original data | |
69 | theta = angle(-sum(rNoCompSampled .^ M)) / M; | |
70 | %% if theta approx +pi/M, wrap to -pi/M | |
71 | if abs(theta - pi / M) / (pi / M) < 0.1 | |
72 | theta = -pi / M; | |
73 | end | |
74 | rNoCompSampled = rNoCompSampled .* exp(-j * theta); | |
75 | ||
76 | ||
77 | %% Not entirely sure why, but after using FFT instead of time-domain | |
78 | %% convolution for simulating CD, we now need to do the same rotation | |
79 | %% for rSampled as well, but this time with a positive rotation. | |
80 | theta = angle(-sum(rSampled .^ M)) / M; | |
81 | if abs(theta + pi / M) / (pi / M) < 0.1 | |
82 | theta = +pi / M; | |
83 | end | |
84 | rSampled = rSampled .* exp(-1j * theta); | |
85 | ||
86 | %% adaptive filter | |
87 | adaptFilterOut = adaptiveCMA(rSampled); | |
88 | ||
89 | demodAdapt = pskdemod(adaptFilterOut, M, pi / M, 'gray'); | |
90 | [~, ber(i)] = biterr(data, demodAdapt) | |
91 | end | |
92 | ||
93 | return | |
94 | ||
95 | ||
96 | figure(1); | |
97 | clf; | |
98 | ||
99 | %% Plot simulated results | |
100 | semilogy(power_dBm, ber, 'Color', [0, 0.6, 0], 'LineWidth', 2); | |
101 | hold on; | |
102 | ||
103 | title({'CD + Kerr + CD compensation', ... | |
104 | strcat(['$D = 17$ ps/(nm km), $z = ', num2str(z), '$ km'])}); | |
105 | grid on; | |
106 | %%xlabel('$E_b/N_0$ (dB)'); | |
107 | xlabel('Optical power (dBm)'); | |
108 | ylabel('BER'); | |
109 | ||
110 | formatFigure; |