| 1 | numSymbs = 5e5; |
| 2 | M = 4; |
| 3 | |
| 4 | Rsym = 2.5e10; % symbol rate (sym/sec) |
| 5 | |
| 6 | rolloff = 0.25; |
| 7 | span = 6; % filter span |
| 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 | 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); |
| 28 | berNoComp = zeros(1, plotlen); |
| 29 | berAdapt = zeros(1, plotlen); |
| 30 | berMatlabAdapt = zeros(1, plotlen); |
| 31 | |
| 32 | data = randi([0 M - 1], numSymbs, 1); |
| 33 | modData = pskmod(data, M, pi / M, 'gray'); |
| 34 | x = txFilter(modData, rolloff, span, sps); |
| 35 | |
| 36 | %% Simulate chromatic dispersion |
| 37 | D = 17; % ps / (nm km) |
| 38 | lambda = 1550; % nm |
| 39 | z = 60;%000; % km |
| 40 | |
| 41 | usingFFT = 1 |
| 42 | xCD = chromaticDispersion_FFT(x, D, lambda, z, Tsamp); |
| 43 | %%xCD = normalizeEnergy(xCD, numSymbs, 1); |
| 44 | %%xCD = x; |
| 45 | |
| 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'); |
| 51 | %%y = xCD; |
| 52 | |
| 53 | r = rxFilter(y, rolloff, span, sps); |
| 54 | rCDComp = CDCompensation(r, D, lambda, z, Tsamp); |
| 55 | rCDComp = normalizeEnergy(rCDComp, numSymbs*sps, 1); |
| 56 | |
| 57 | rSampled = rCDComp(sps*span/2+1:sps:(numSymbs+span/2)*sps); |
| 58 | rNoCompSampled = r(sps*span/2+1:sps:(numSymbs+span/2)*sps); |
| 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 |
| 66 | rNoCompSampled = rNoCompSampled .* exp(-j * theta); |
| 67 | |
| 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 | |
| 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'); |
| 87 | |
| 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); |
| 92 | |
| 93 | %{ |
| 94 | if EbN0_db(i) == 14 |
| 95 | figure(1); |
| 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); |
| 102 | %%scatterplot(modData); |
| 103 | %%title('Original constellation'); |
| 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'); |
| 112 | %scatterplot(matlabEq); |
| 113 | %title('Matlab equalizer'); |
| 114 | ber(i) |
| 115 | %berNoComp(i) |
| 116 | %berAdapt(i) |
| 117 | %berMatlabAdapt(i) |
| 118 | end |
| 119 | %} |
| 120 | end |
| 121 | |
| 122 | figure(1); |
| 123 | clf; |
| 124 | |
| 125 | %% Plot simulated results |
| 126 | semilogy(EbN0_db, ber, 'r', 'LineWidth', 2); |
| 127 | hold on; |
| 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); |
| 131 | |
| 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', ... |
| 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'])}); |
| 144 | grid on; |
| 145 | xlabel('$E_b/N_0$ (dB)'); |
| 146 | ylabel('BER'); |
| 147 | |
| 148 | formatFigure; |