CD_AWGN.m

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