kerr.m

   1 numSymbs = 2^16;
   2 M = 4;
   3
   4 Rsym = 2.5e10; % symbol rate (sym/sec)
   5 Tsym = 1 / Rsym; % symbol period (sec)
   6
   7 rolloff = 0.25;
   8 span = 6; % filter span
   9 sps = 8; % samples per symbol
  10
  11 fs = Rsym * sps; % sampling freq (Hz)
  12 Tsamp = 1 / fs;
  13
  14 t = (0 : 1 / fs : numSymbs / Rsym + (1.5 * span * sps - 1) / fs).';
  15
  16 power_dBm = -6:1:4;
  17 %%power_dBm = 0;
  18 power = 10 .^ (power_dBm / 10) * 1e-3; % watts
  19
  20 Es = power * Tsym; % joules
  21 Eb = Es / log2(M); % joules
  22
  23 N0ref_db = 10; % Eb/N0 at power = 1mW
  24 %% Fix N0, such that Eb/N0 = N0ref_db at power = 1mW
  25 N0 = 1e-3 * Tsym / (log2(M) * 10 ^ (N0ref_db / 10)); % joules
  26 %% At current settings, N0 = 0.002 pJ
  27
  28 plotlen = length(power);
  29
  30 ber = zeros(1, plotlen);
  31
  32 data = randi([0 M - 1], numSymbs, 1);
  33 modData = dpskmod(data, M, 0, 'gray');
  34 %%modData = pskmod(data, M, pi/4, 'gray');
  35
  36
  37 %% Chromatic dispersion
  38 D = 17; % ps / (nm km)
  39 lambda = 1550; % nm
  40 z = 100; % km
  41
  42
  43 TsampOrig = Tsamp;
  44
  45 for i = 1:plotlen
  46   sps = 8;
  47   Tsamp = TsampOrig;
  48
  49   snr = Es(i) / sps / N0;
  50   snr_dB = 10 * log10(snr);
  51
  52   x = txFilter(modData, rolloff, span, sps);
  53   %% Now, sum(abs(x) .^ 2) / length(x) should be 1.
  54   %% We can set its power simply by multiplying.
  55   x = sqrt(power(i)) * x;
  56
  57   %% We can now do split-step Fourier.
  58   gamma = 1.2; % watt^-1 / km
  59
  60
  61   xCDKerr = splitstepfourier(x, D, lambda, z, Tsamp, gamma);
  62
  63   y = awgn(xCDKerr, snr_dB, 'measured', 'db');
  64   %y = xCDKerr;
  65
  66   r = rxFilter(y, rolloff, span, sps);
  67   sps = 2;
  68   Tsamp = Tsamp * 4;
  69
  70   rCDComp = CDCompensation(r, D, lambda, z, Tsamp);
  71   rCDComp = normalizeEnergy(rCDComp, numSymbs * sps, 1);
  72
  73   rSampled = rCDComp(2:2:end);
  74
  75   %% adaptive filter
  76   [adaptFilterOut, convergeIdx] = adaptiveCMA(rSampled);
  77
  78   demod = dpskdemod(adaptFilterOut, M, 0, 'gray');
  79   %%demod = pskdemod(adaptFilterOut, M, pi/4, 'gray');
  80
  81   if convergeIdx < Inf
  82     [~, ber(i)] = biterr(data(convergeIdx:end), demod(convergeIdx:end));
  83   else
  84     [~, ber(i)] = biterr...
  85                     (data(ceil(0.8*numSymbs):end), ...
  86                      demod(ceil(0.8*numSymbs):end));
  87   end
  88 end
  89
  90 ber
  91
  92
  93 figure;
  94 clf;
  95
  96 %% Plot simulated results
  97 qp = 20 * log10(erfcinv(2*ber)*sqrt(2));
  98 plot(power_dBm, qp, 'Color', [0, 0.6, 0], 'LineWidth', 2);
  99 hold on;
 100
 101 title({'CD + Kerr + CD compensation', ...
 102        strcat(['$D = 17$ ps/(nm km), $z = ', num2str(z), '$ km'])});
 103 grid on;
 104 xlabel('Optical power (dBm)');
 105 ylabel('$20 \log_{10}\left(\sqrt{2}\mathrm{erfc}^{-1}(2 BER)\right)$');
 106
 107 formatFigure;