| 1 | M = 4; |
| 2 | numSymbs = 5e5; |
| 3 | |
| 4 | Rsym = 2.5e10; % symbol rate (sym/sec) |
| 5 | |
| 6 | span = 6; % Tx/Rx filter span |
| 7 | rolloff = 0.25; % Tx/Rx RRC rolloff |
| 8 | sps = 2; % samples per symbol |
| 9 | |
| 10 | fs = Rsym * sps; % sampling freq (Hz) |
| 11 | Tsamp = 1 / fs; |
| 12 | |
| 13 | t = (0 : 1 / fs : numSymbs / Rsym + (1.5 * span * sps - 1) / fs).'; |
| 14 | |
| 15 | data = randi([0 M - 1], numSymbs, 1); |
| 16 | modData = pskmod(data, M, pi / M, 'gray'); |
| 17 | x = txFilter(modData, rolloff, span, sps); |
| 18 | |
| 19 | x = normalizeEnergy(x, numSymbs*sps, 1); |
| 20 | |
| 21 | %% Simulate chromatic dispersion |
| 22 | D = 17; % ps / (nm km) |
| 23 | lambda = 1550; % nm |
| 24 | z = 5000 % km |
| 25 | |
| 26 | [xCD, xCDkstart] = chromaticDispersion_FFT(x, D, lambda, z, Tsamp); |
| 27 | |
| 28 | EbN0_db = 8; |
| 29 | snr = EbN0_db + 10 * log10(log2(M)) - 10 * log10(sps); |
| 30 | |
| 31 | %%y = awgn(xCD, snr, 'measured'); |
| 32 | y = xCD; |
| 33 | |
| 34 | r = rxFilter(y, rolloff, span, sps); |
| 35 | [rCDComp, CDCompkstart] = CDCompensation(r, D, lambda, z, Tsamp); |
| 36 | rCDComp = normalizeEnergy(rCDComp, numSymbs*sps, 1); |
| 37 | |
| 38 | rSampled = rCDComp(sps*span/2+1:sps:(numSymbs+span/2)*sps); |
| 39 | rNoCompSa = r(sps*span/2+1:sps:(numSymbs+span/2)*sps); |
| 40 | |
| 41 | %% if no CD comp, then rotate constellation. Use: |
| 42 | theta = angle(-sum(rNoCompSa .^ M)) / M; |
| 43 | %% if theta approx +pi/M, wrap to -pi/M |
| 44 | if abs(theta - pi / M) / (pi / M) < 0.1 |
| 45 | theta = -pi / M; |
| 46 | end |
| 47 | rNoCompSa = rNoCompSa .* exp(-j * theta); |
| 48 | |
| 49 | |
| 50 | %% Not entirely sure why, but after using FFT instead of time-domain |
| 51 | %% convolution for simulating CD, we now need to do the same rotation |
| 52 | %% for rSampled as well, but this time with a positive rotation. |
| 53 | theta = angle(-sum(rSampled .^ M)) / M; |
| 54 | if abs(theta + pi / M) / (pi / M) < 0.1 |
| 55 | theta = +pi / M; |
| 56 | end |
| 57 | rSampled = rSampled .* exp(-1j * theta); |
| 58 | |
| 59 | |
| 60 | |
| 61 | %%rAdaptEq = adaptiveCMA(rSampled); |
| 62 | %{ |
| 63 | %% Compare original signal and compensated signal |
| 64 | figure(101); |
| 65 | clf; |
| 66 | tsym = t(sps*span/2+1:sps:(numSymbs+span/2)*sps); |
| 67 | subplot(211); |
| 68 | plot(t(1:length(x)), real(normalizeEnergy(x, numSymbs*sps, 1)), 'b'); |
| 69 | hold on |
| 70 | plot(t(1:length(x)), real(normalizeEnergy(yCDComp(1:length(x)), numSymbs*sps, 1)), 'r'); |
| 71 | plot(tsym, real(rAdaptEq), 'x', 'Color', [0, 0.6, 0], 'LineWidth', 2); |
| 72 | hold off; |
| 73 | title('Real part'); |
| 74 | legend('original', 'dispersion compensated', 'CMA equalized samples'); |
| 75 | axis([t(6000*sps+1) t(6000*sps+150) -Inf +Inf]); |
| 76 | subplot(212); |
| 77 | plot(t(1:length(x)), imag(normalizeEnergy(x, numSymbs*sps, 1)), 'b'); |
| 78 | hold on; |
| 79 | plot(t(1:length(x)), imag(normalizeEnergy(yCDComp(1:length(x)), numSymbs*sps, 1)), 'r'); |
| 80 | plot(tsym, imag(rAdaptEq), 'x', 'Color', [0, 0.6, 0], 'LineWidth', 2); |
| 81 | hold off; |
| 82 | title('Imag part'); |
| 83 | axis([t(6000*sps+1) t(6000*sps+150) -Inf +Inf]); |
| 84 | |
| 85 | scatterplot(modData); |
| 86 | formatFigure; |
| 87 | %title('Constellation of original modulation', 'interpreter', 'latex'); |
| 88 | xlabel('In-Phase', 'interpreter', 'latex'); |
| 89 | %scatterplot(rSampled); |
| 90 | %title('Constellation of matched filter output'); |
| 91 | scatterplot(rNoCompSa); |
| 92 | title('Constellation of dispersed signal', 'interpreter', 'latex'); |
| 93 | scatterplot(rAdaptEq); |
| 94 | title('Constellation of adaptive filter output'); |
| 95 | %} |
| 96 | demodData = pskdemod(rSampled, M, pi / M, 'gray'); |
| 97 | %%demodAdapt = pskdemod(rAdaptEq, M, pi / M, 'gray'); |
| 98 | |
| 99 | [~, ber] = biterr(data, demodData) |
| 100 | %[~, berNoComp] = biterr(data, pskdemod(rNoCompSa, M, pi/M, 'gray')) |
| 101 | %[~, ber] = biterr(data, demodAdapt) |