4 Rsym = 2.5e10; % symbol rate (sym/sec)
7 span = 6; % filter span
8 sps = 2; % samples per symbol
10 fs = Rsym * sps; % sampling freq (Hz)
13 t = (0 : 1 / fs : numSymbs / Rsym + (1.5 * span * sps - 1) / fs).';
16 EbN0 = 10 .^ (EbN0_db ./ 10);
22 EsN0 = EbN0 .* log2(M);
23 EsN0_db = 10 .* log10(EsN0);
25 plotlen = length(EbN0);
27 ber = zeros(1, plotlen);
28 berNoComp = zeros(1, plotlen);
29 berAdapt = zeros(1, plotlen);
30 berMatlabAdapt = zeros(1, plotlen);
32 data = randi([0 M - 1], numSymbs, 1);
33 modData = pskmod(data, M, pi / M, 'gray');
34 x = txFilter(modData, rolloff, span, sps);
36 %% Simulate chromatic dispersion
37 D = 17; % ps / (nm km)
42 xCD = chromaticDispersion_FFT(x, D, lambda, z, Tsamp);
43 %%xCD = normalizeEnergy(xCD, numSymbs, 1);
47 snr = EbN0_db(i) + 10 * log10(log2(M)) - 10 * log10(sps);
48 noiseEnergy = 10 ^ (-snr / 10);
50 y = awgn(xCD, snr, 'measured');
53 r = rxFilter(y, rolloff, span, sps);
54 rCDComp = CDCompensation(r, D, lambda, z, Tsamp);
55 rCDComp = normalizeEnergy(rCDComp, numSymbs*sps, 1);
57 rSampled = rCDComp(sps*span/2+1:sps:(numSymbs+span/2)*sps);
58 rNoCompSampled = r(sps*span/2+1:sps:(numSymbs+span/2)*sps);
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
66 rNoCompSampled = rNoCompSampled .* exp(-j * theta);
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
76 rSampled = rSampled .* exp(-1j * theta);
81 adaptFilterOut = adaptiveCMA(rSampled);
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');
88 [bitErrors, ber(i)] = biterr(data, demodData);
89 [bitErrors, berNoComp(i)] = biterr(data, demodNoComp);
90 [~, berAdapt(i)] = biterr(data, demodAdapt);
91 %%[~, berMatlabAdapt(i)] = biterr(data, demodMatlabAdapt);
96 scatterplot(normalizeEnergy(rSampled, numSymbs, 1));
98 title('Constellation after CD comp.', 'interpreter', 'latex');
99 xlabel('In-Phase', 'interpreter', 'latex');
100 ylabel('Quadrature', 'interpreter', 'latex');
101 set(gca, 'FontSize', 18);
102 %%scatterplot(modData);
103 %%title('Original constellation');
104 scatterplot(normalizeEnergy(rNoCompSampled, numSymbs, 1));
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');
112 %scatterplot(matlabEq);
113 %title('Matlab equalizer');
125 %% Plot simulated results
126 semilogy(EbN0_db, ber, 'r', 'LineWidth', 2);
128 semilogy(EbN0_db, berNoComp, 'm', 'LineWidth', 2);
129 semilogy(EbN0_db, berAdapt, 'Color', [0, 0.6, 0], 'LineWidth', 2);
130 %%%semilogy(EbN0_db, berMatlabAdapt, 'c', 'LineWidth', 1.4);
132 theoreticalPSK(EbN0_db, M, 'b', 'LineWidth', 1);
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', ...
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'])});
145 xlabel('$E_b/N_0$ (dB)');