4 Rsym = 2.5e10; % symbol rate (sym/sec)
7 span = 6; % filter span
8 sps = 4; % 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)
41 xCD = chromaticDispersion(x, D, lambda, z, Tsamp);
42 xCD = normalizeEnergy(xCD, numSymbs, 1);
45 snr = EbN0_db(i) + 10 * log10(log2(M)) - 10 * log10(sps);
46 noiseEnergy = 10 ^ (-snr / 10);
48 y = awgn(xCD, snr, 'measured');
50 yCDComp = CDCompensation(y, D, lambda, z, Tsamp);
52 r = rxFilter(yCDComp, rolloff, span, sps);
53 rNoComp = rxFilter(y, rolloff, span, sps);
55 %r = normalizeEnergy(r, numSymbs, 1 + noiseEnergy);
57 rSampled = r(sps*span/2+1:sps:(numSymbs + span/2) * sps);
58 rNoCompSampled = rNoComp(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
67 rNoCompSampled = rNoCompSampled .* exp(-j * theta);
70 adaptFilterOut = adaptiveCMA(rSampled);
72 demodData = pskdemod(rSampled, M, pi / M, 'gray');
73 demodNoComp = pskdemod(rNoCompSampled, M, pi / M, 'gray');
74 demodAdapt = pskdemod(adaptFilterOut, M, pi / M, 'gray');
75 %%demodMatlabAdapt = pskdemod(matlabEq, M, pi / M, 'gray');
77 [bitErrors, ber(i)] = biterr(data, demodData);
78 [bitErrors, berNoComp(i)] = biterr(data, demodNoComp);
79 [~, berAdapt(i)] = biterr(data, demodAdapt);
80 %%[~, berMatlabAdapt(i)] = biterr(data, demodMatlabAdapt);
85 scatterplot(rSampled);
86 title('Constellation after CD compensation');
87 %%scatterplot(modData);
88 %%title('Original constellation');
89 scatterplot(rNoCompSampled);
90 title('Constellation without CD compensation');
91 scatterplot(adaptFilterOut);
92 title('Constellation with CD compensation and adaptive filter');
93 %scatterplot(matlabEq);
94 %title('Matlab equalizer');
106 %% Plot simulated results
107 semilogy(EbN0_db, ber, 'r', 'LineWidth', 2);
109 %%semilogy(EbN0_db, berNoComp, 'g', 'LineWidth', 2);
110 semilogy(EbN0_db, berAdapt, 'm', 'LineWidth', 1.4);
111 %%%semilogy(EbN0_db, berMatlabAdapt, 'c', 'LineWidth', 1.4);
113 theoreticalPSK(EbN0_db, M, 'b', 'LineWidth', 1);
114 legend({'CD + AWGN + CD comp.', 'CD + AWGN + CD comp.~+ CMA', ...
115 'Theoretical AWGN'}, 'Location', 'southwest');
117 title(strcat(num2str(M), '-PSK with chromatic dispersion and compensation'));
119 xlabel('$E_b/N_0$ (dB)');