kerr1Signal.m

   1 numSymbs = 5e5;
   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 = 2; % 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
  17 %%power_dBm = -3:0.2:4;
  18 power_dBm = [0];
  19 power = 10 .^ (power_dBm / 10) * 1e-3; % watts
  20
  21 Es = power * Tsym; % joules
  22 Eb = Es / log2(M); % joules
  23
  24 N0ref_db = 10; % Eb/N0 at power = 1mW
  25 %% Fix N0, such that Eb/N0 = N0ref_db at power = 1mW
  26 N0 = 1e-3 * Tsym / (log2(M) * 10 ^ (N0ref_db / 10)); % joules
  27
  28
  29 plotlen = length(power);
  30
  31 ber = zeros(1, plotlen);
  32
  33 data = randi([0 M - 1], numSymbs, 1);
  34 modData = pskmod(data, M, pi / M, 'gray');
  35
  36
  37 %% Chromatic dispersion
  38 D = 17; % ps / (nm km)
  39 lambda = 1550; % nm
  40 z = 600; % km
  41
  42
  43 for i = 1:plotlen
  44   snr = Es(i) / sps / N0;
  45   snr_dB = 10 * log10(snr);
  46
  47   x = txFilter(modData, rolloff, span, sps);
  48   %% Now, sum(abs(x) .^ 2) / length(x) should be 1.
  49   %% We can set its power simply by multiplying.
  50   x = sqrt(power(i)) * x;
  51
  52   %% We can now do split-step Fourier.
  53   gamma = 1.2; % watt^-1 / km
  54   %%stepnum = round(40 * z * gamma); % Nonlinear Fiber optics, App B
  55   stepnum = 100;
  56   xCD = splitstepfourier(x, D, lambda, z, Tsamp, gamma, stepnum);
  57
  58   y = awgn(xCD, snr, power(i), 'linear');
  59   %%y = xCD;
  60
  61   r = rxFilter(y, rolloff, span, sps);
  62   rCDComp = CDCompensation(r, D, lambda, z, Tsamp);
  63   rCDComp = normalizeEnergy(rCDComp, numSymbs*sps, 1);
  64
  65   rSampled = rCDComp(sps*span/2+1:sps:(numSymbs+span/2)*sps);
  66   rNoCompSampled = r(sps*span/2+1:sps:(numSymbs+span/2)*sps);
  67
  68   %% rotate rNoCompSampled to match original data
  69   theta = angle(-sum(rNoCompSampled .^ M)) / M;
  70   %% if theta approx +pi/M, wrap to -pi/M
  71   if abs(theta - pi / M) / (pi / M) < 0.1
  72     theta = -pi / M;
  73   end
  74   rNoCompSampled = rNoCompSampled .* exp(-j * theta);
  75
  76
  77   %% Not entirely sure why, but after using FFT instead of time-domain
  78   %% convolution for simulating CD, we now need to do the same rotation
  79   %% for rSampled as well, but this time with a positive rotation.
  80   theta = angle(-sum(rSampled .^ M)) / M;
  81   if abs(theta + pi / M) / (pi / M) < 0.1
  82     theta = +pi / M;
  83   end
  84   rSampled = rSampled .* exp(-1j * theta);
  85
  86   %% adaptive filter
  87   adaptFilterOut = adaptiveCMA(rSampled);
  88
  89   demodAdapt = pskdemod(adaptFilterOut, M, pi / M, 'gray');
  90   [~, ber(i)] = biterr(data, demodAdapt)
  91 end
  92
  93 return
  94
  95
  96 figure(1);
  97 clf;
  98
  99 %% Plot simulated results
 100 semilogy(power_dBm, ber, 'Color', [0, 0.6, 0], 'LineWidth', 2);
 101 hold on;
 102
 103 title({'CD + Kerr + CD compensation', ...
 104        strcat(['$D = 17$ ps/(nm km), $z = ', num2str(z), '$ km'])});
 105 grid on;
 106 %%xlabel('$E_b/N_0$ (dB)');
 107 xlabel('Optical power (dBm)');
 108 ylabel('BER');
 109
 110 formatFigure;