channel.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, 0, 'gray');
  35 for i = 2:numSymbs
  36   modData(i) = modData(i) * modData(i-1);
  37 end
  38
  39
  40 %% Chromatic dispersion
  41 D = 17; % ps / (nm km)
  42 lambda = 1550; % nm
  43 z = 100; % km
  44
  45
  46 linewidthTx = 0; % Hz
  47 linewidthLO = 1e6; % Hz
  48
  49
  50 TsampOrig = Tsamp;
  51
  52 x_P1 = txFilter(modData, rolloff, span, sps);
  53
  54
  55 for i = 1:plotlen
  56   sps = 8;
  57   Tsamp = TsampOrig;
  58
  59   snr = Es(i) / sps / N0;
  60   snr_dB = 10 * log10(snr);
  61
  62   %%x = txFilter(modData, rolloff, span, sps);
  63   %% Now, sum(abs(x) .^ 2) / length(x) should be 1.
  64   %% We can set its power simply by multiplying.
  65   x = sqrt(power(i)) * x_P1;
  66
  67   %% We can now do split-step Fourier.
  68   gamma = 1.2; % watt^-1 / km
  69
  70
  71   xCDKerr = splitstepfourier(x, D, lambda, z, Tsamp, gamma);
  72
  73   xpn = phaseNoise(xCDKerr, linewidthTx, linewidthLO, Tsamp);
  74
  75   y = awgn(xpn, snr_dB, 'measured', 'db');
  76   %y = xCDKerr;
  77
  78   r = rxFilter(y, rolloff, span, sps);
  79   sps = 2;
  80   Tsamp = Tsamp * 4;
  81
  82   rCDComp = CDCompensation(r, D, lambda, z, Tsamp);
  83   rCDComp = normalizeEnergy(rCDComp, numSymbs * sps, 1);
  84
  85   rSampled = rCDComp(2:2:end);
  86
  87   %% adaptive filter
  88   [adaptFilterOut, convergeIdx] = adaptiveCMA(rSampled);
  89
  90   pncorr = phaseNoiseCorr(adaptFilterOut, M, 0, 40).';
  91
  92   demodAdapt = pskdemod(pncorr, M, 0, 'gray');
  93   remod = pskmod(demodAdapt, M, 0, 'gray');
  94   delayed = [1; remod(1:end-1)];
  95   demod = pskdemod(remod .* conj(delayed), M, 0, 'gray');
  96
  97   if convergeIdx < Inf
  98     [~, ber(i)] = biterr(data(convergeIdx:end), demod(convergeIdx:end));
  99   else
 100     [~, ber(i)] = biterr...
 101                     (data(ceil(0.8*numSymbs):end), ...
 102                      demod(ceil(0.8*numSymbs):end));
 103   end
 104 end
 105
 106 ber
 107
 108
 109 figure(1);
 110 clf;
 111
 112 %% Plot simulated results
 113 qp = 20 * log10(erfcinv(2*ber)*sqrt(2));
 114 plot(power_dBm, qp, 'Color', [0, 0.6, 0], 'LineWidth', 2);
 115 hold on;
 116
 117 title({'CD + Kerr + CD compensation', ...
 118        strcat(['$D = 17$ ps/(nm km), $z = ', num2str(z), '$ km'])});
 119 grid on;
 120 xlabel('Optical power (dBm)');
 121 ylabel('$20 \log_{10}\left(\sqrt{2}\mathrm{erfc}^{-1}(2 BER)\right)$');
 122
 123 formatFigure;