Commit | Line | Data |
---|---|---|
f9a73e9e | 1 | numSymbs = 5e5; |
1eeb62fb AIL |
2 | M = 4; |
3 | ||
4 | Rsym = 2.5e10; % symbol rate (sym/sec) | |
5 | ||
6 | rolloff = 0.25; | |
7 | span = 6; % filter span | |
f9a73e9e | 8 | sps = 2; % samples per symbol |
1eeb62fb AIL |
9 | |
10 | fs = Rsym * sps; % sampling freq (Hz) | |
11 | Tsamp = 1 / fs; | |
12 | ||
5e9be3c4 | 13 | t = (0 : 1 / fs : numSymbs / Rsym + (1.5 * span * sps - 1) / fs).'; |
1eeb62fb AIL |
14 | |
15 | EbN0_db = 0:0.2:14; | |
16 | EbN0 = 10 .^ (EbN0_db ./ 10); | |
17 | ||
18 | Es = 1; | |
19 | Eb = Es / log2(M); | |
20 | N0 = Eb ./ EbN0; | |
21 | ||
22 | EsN0 = EbN0 .* log2(M); | |
23 | EsN0_db = 10 .* log10(EsN0); | |
24 | ||
25 | plotlen = length(EbN0); | |
26 | ||
27 | ber = zeros(1, plotlen); | |
5e9be3c4 AIL |
28 | berNoComp = zeros(1, plotlen); |
29 | berAdapt = zeros(1, plotlen); | |
30 | berMatlabAdapt = zeros(1, plotlen); | |
1eeb62fb AIL |
31 | |
32 | data = randi([0 M - 1], numSymbs, 1); | |
5e9be3c4 | 33 | modData = pskmod(data, M, pi / M, 'gray'); |
1eeb62fb AIL |
34 | x = txFilter(modData, rolloff, span, sps); |
35 | ||
36 | %% Simulate chromatic dispersion | |
5e9be3c4 | 37 | D = 17; % ps / (nm km) |
1eeb62fb | 38 | lambda = 1550; % nm |
f9a73e9e | 39 | z = 60;%000; % km |
1eeb62fb | 40 | |
f9a73e9e AIL |
41 | usingFFT = 1 |
42 | xCD = chromaticDispersion_FFT(x, D, lambda, z, Tsamp); | |
43 | %%xCD = normalizeEnergy(xCD, numSymbs, 1); | |
44 | %%xCD = x; | |
1eeb62fb | 45 | |
1eeb62fb AIL |
46 | for i = 1:plotlen |
47 | snr = EbN0_db(i) + 10 * log10(log2(M)) - 10 * log10(sps); | |
48 | noiseEnergy = 10 ^ (-snr / 10); | |
49 | ||
50 | y = awgn(xCD, snr, 'measured'); | |
f9a73e9e | 51 | %%y = xCD; |
1eeb62fb | 52 | |
f9a73e9e AIL |
53 | r = rxFilter(y, rolloff, span, sps); |
54 | rCDComp = CDCompensation(r, D, lambda, z, Tsamp); | |
55 | rCDComp = normalizeEnergy(rCDComp, numSymbs*sps, 1); | |
1eeb62fb | 56 | |
f9a73e9e AIL |
57 | rSampled = rCDComp(sps*span/2+1:sps:(numSymbs+span/2)*sps); |
58 | rNoCompSampled = r(sps*span/2+1:sps:(numSymbs+span/2)*sps); | |
5e9be3c4 AIL |
59 | |
60 | %% rotate rNoCompSampled to match original data | |
61 | theta = angle(-sum(rNoCompSampled .^ M)) / M; | |
62 | %% if theta approx +pi/M, wrap to -pi/M | |
63 | if abs(theta - pi / M) / (pi / M) < 0.1 | |
64 | theta = -pi / M; | |
65 | end | |
5e9be3c4 AIL |
66 | rNoCompSampled = rNoCompSampled .* exp(-j * theta); |
67 | ||
f9a73e9e AIL |
68 | |
69 | %% Not entirely sure why, but after using FFT instead of time-domain | |
70 | %% convolution for simulating CD, we now need to do the same rotation | |
71 | %% for rSampled as well, but this time with a positive rotation. | |
72 | theta = angle(-sum(rSampled .^ M)) / M; | |
73 | if abs(theta + pi / M) / (pi / M) < 0.1 | |
74 | theta = +pi / M; | |
75 | end | |
76 | rSampled = rSampled .* exp(-1j * theta); | |
77 | ||
78 | ||
79 | ||
5e9be3c4 AIL |
80 | %% adaptive filter |
81 | adaptFilterOut = adaptiveCMA(rSampled); | |
82 | ||
83 | demodData = pskdemod(rSampled, M, pi / M, 'gray'); | |
84 | demodNoComp = pskdemod(rNoCompSampled, M, pi / M, 'gray'); | |
85 | demodAdapt = pskdemod(adaptFilterOut, M, pi / M, 'gray'); | |
86 | %%demodMatlabAdapt = pskdemod(matlabEq, M, pi / M, 'gray'); | |
1eeb62fb AIL |
87 | |
88 | [bitErrors, ber(i)] = biterr(data, demodData); | |
5e9be3c4 AIL |
89 | [bitErrors, berNoComp(i)] = biterr(data, demodNoComp); |
90 | [~, berAdapt(i)] = biterr(data, demodAdapt); | |
91 | %%[~, berMatlabAdapt(i)] = biterr(data, demodMatlabAdapt); | |
92 | ||
f9a73e9e AIL |
93 | %{ |
94 | if EbN0_db(i) == 14 | |
5e9be3c4 | 95 | figure(1); |
f9a73e9e AIL |
96 | scatterplot(normalizeEnergy(rSampled, numSymbs, 1)); |
97 | formatFigure; | |
98 | title('Constellation after CD comp.', 'interpreter', 'latex'); | |
99 | xlabel('In-Phase', 'interpreter', 'latex'); | |
100 | ylabel('Quadrature', 'interpreter', 'latex'); | |
101 | set(gca, 'FontSize', 18); | |
5e9be3c4 AIL |
102 | %%scatterplot(modData); |
103 | %%title('Original constellation'); | |
f9a73e9e AIL |
104 | scatterplot(normalizeEnergy(rNoCompSampled, numSymbs, 1)); |
105 | formatFigure; | |
106 | title('Constellation without CD comp.', 'interpreter', 'latex'); | |
107 | xlabel('In-Phase', 'interpreter', 'latex'); | |
108 | ylabel('Quadrature', 'interpreter', 'latex'); | |
109 | set(gca, 'FontSize', 18); | |
110 | %scatterplot(adaptFilterOut); | |
111 | %title('Constellation with CD compensation and adaptive filter'); | |
5e9be3c4 AIL |
112 | %scatterplot(matlabEq); |
113 | %title('Matlab equalizer'); | |
114 | ber(i) | |
115 | %berNoComp(i) | |
f9a73e9e AIL |
116 | %berAdapt(i) |
117 | %berMatlabAdapt(i) | |
5e9be3c4 | 118 | end |
f9a73e9e | 119 | %} |
1eeb62fb AIL |
120 | end |
121 | ||
122 | figure(1); | |
123 | clf; | |
124 | ||
125 | %% Plot simulated results | |
126 | semilogy(EbN0_db, ber, 'r', 'LineWidth', 2); | |
127 | hold on; | |
f9a73e9e AIL |
128 | semilogy(EbN0_db, berNoComp, 'm', 'LineWidth', 2); |
129 | semilogy(EbN0_db, berAdapt, 'Color', [0, 0.6, 0], 'LineWidth', 2); | |
5e9be3c4 | 130 | %%%semilogy(EbN0_db, berMatlabAdapt, 'c', 'LineWidth', 1.4); |
1eeb62fb AIL |
131 | |
132 | theoreticalPSK(EbN0_db, M, 'b', 'LineWidth', 1); | |
f9a73e9e AIL |
133 | %%legend({'CD + AWGN + CD comp.', 'CD + AWGN + CD comp.~+ CMA', ... |
134 | %% 'Theoretical AWGN'}, 'Location', 'southwest'); | |
135 | %%legend({'CD + AWGN + CD comp.', 'CD + AWGN', 'Theoretical AWGN'}, ... | |
136 | %% 'Location', 'southwest'); | |
137 | legend({'CD + AWGN + CD comp.', 'CD + AWGN', ... | |
138 | 'CD + AWGN + CD comp.~+ CMA', 'Theoretical AWGN'}, 'Location', ... | |
139 | 'Southwest'); | |
140 | ||
141 | %%title(strcat(num2str(M), '-PSK with chromatic dispersion and compensation')); | |
142 | title({'QPSK with chromatic dispersion and compensation', ... | |
143 | strcat(['$D = 17$ ps/(nm km), $z = ', num2str(z), '$ km'])}); | |
1eeb62fb AIL |
144 | grid on; |
145 | xlabel('$E_b/N_0$ (dB)'); | |
146 | ylabel('BER'); | |
147 | ||
148 | formatFigure; |