main.m

   1 numSymbs = 2^16;
   2 M = 4; % QPSK
   3
   4 Rsym = 28e9; % symbol rate (sym/sec)
   5
   6 %% zs: array of distances z to be simulated
   7 % Example: zs = 42; zs = 40:10:100; zs = [300, 500, 1000];
   8 zs = 42;
   9 plotlen = length(zs);
  10
  11
  12 %% Tx RRC filter properties
  13 rolloff = 0.25;
  14 span = 6; % filter span
  15 sps = 16; % samples per symbol
  16
  17 %% Sampling frequency
  18 fs = Rsym * sps; % Hz
  19 Tsamp = 1 / fs; % s
  20 % t: time vector, s
  21 t = (0 : 1 / fs : numSymbs / Rsym - 1 / fs).';
  22
  23
  24 %% Chromatic dispersion
  25 D = 17; % ps / (nm km)
  26 lambda = 1550; % nm
  27
  28 %% Laser phase noise
  29 linewidthTx = 0; % Hz
  30 linewidthLO = 1e6; % Hz
  31
  32 %% Kerr effect / SSFS parameters
  33 gamma = 1.2; % watt^-1 / km
  34 alpha = 0.2; % dB/km
  35 dz = 2; % Step size, km
  36
  37 %% Polarization state rotation parameters
  38 rot_omega = 1e3; % rad/s
  39 rot_phi = 2; % rad
  40
  41 %% Launch power, per wavelength channel
  42 power_dBm = 0;
  43 power = 10 .^ (power_dBm / 10) * 1e-3; % watts
  44
  45 %% WDM properties
  46 wavelength_channels = 3;
  47 dw = 2 * pi * 50e9; % channel spacing (rad/s)
  48
  49 %% Shot noise
  50 hc = 6.62607015e-34 * 299792458; % J m
  51 Eperphoton = hc / (lambda * 1e-9); % J
  52
  53
  54 %% Stores result to be plotted
  55 ber = zeros(plotlen, 1);
  56 if plotlen > 1
  57   fig = figure; hold on;
  58 end
  59
  60 %% sps and Tsamp change at Tx/Rx, save these for later.
  61 spsOrig = sps;
  62 TsampOrig = Tsamp;
  63
  64 %% Generate random data for both polarizations
  65 data_x = randi([0, M - 1], numSymbs, wavelength_channels, 'uint8');
  66 data_y = randi([0, M - 1], numSymbs, wavelength_channels, 'uint8');
  67
  68 %% DE-QPSK modulation
  69 modData_x = deqpskmod(data_x);
  70 modData_y = deqpskmod(data_y);
  71
  72 %% Construct waveforms for each channel separately
  73 A_x_wdm = zeros(numSymbs * sps, wavelength_channels);
  74 A_y_wdm = zeros(numSymbs * sps, wavelength_channels);
  75 carriers = zeros(numSymbs * sps, wavelength_channels);
  76
  77 for w = 1 : wavelength_channels
  78   %% Compute frequency offsets:
  79   %                     ___     ___     ___     ___     ___
  80   %       Spectrum     |   |   |   |   |   |   |   |   |   |
  81   %                    |   |   |   |   |   |   |   |   |   |
  82   %                ____|   |___|   |___|   |___|   |___|   |____  --> freq
  83   % channel #            5       3       1       2       4
  84   % ang freq offset    -2dw     -dw      0      +dw    +2dw
  85
  86   if mod(w, 2) == 0
  87     ndw = w / 2 * dw;
  88   else
  89     ndw = (1-w) / 2 * dw;
  90   end
  91   carriers(:, w) = exp(1j * ndw * t);
  92   A_x_wdm(:, w) = txFilter(modData_x(:, w), rolloff, span, sps);
  93   A_y_wdm(:, w) = txFilter(modData_y(:, w), rolloff, span, sps);
  94 end
  95
  96 %% Sum the WDM waveforms with their frequency offsets
  97 A_x = sum(A_x_wdm .* carriers, 2);
  98 A_y = sum(A_y_wdm .* carriers, 2);
  99
 100 %% Clear variables no longer needed to reduce memory usage
 101 clear modData_x modData_y A_x_wdm A_y_wdm;
 102
 103 %% Set launch power. Divide by 2 because half power for each polarization.
 104 A_x = sqrt(power / 2) * A_x;
 105 A_y = sqrt(power / 2) * A_y;
 106
 107 %% Rotate polarization states
 108 A_x = A_x .* cos(rot_omega * t) + ...
 109       A_y .* sin(rot_omega * t) * exp(-1j * rot_phi);
 110 A_y = A_x .* -sin(rot_omega * t) * exp(1j * rot_phi) + ...
 111       A_y .* cos(rot_omega * t);
 112
 113 %% Now loop through each z
 114 for i = 1 : plotlen
 115   z = zs(i);
 116
 117   sps = spsOrig;
 118   Tsamp = TsampOrig;
 119
 120   %% Split-step Fourier
 121   [A_x, A_y] = ssfs(A_x, A_y, D, lambda, z, dz, Tsamp, gamma, alpha);
 122
 123   %% Phase noise
 124   A_x = phaseNoise(A_x, linewidthTx, linewidthLO, Tsamp);
 125   A_y = phaseNoise(A_y, linewidthTx, linewidthLO, Tsamp);
 126
 127   %% Here, only receive the central channel 1.
 128   % For channel n: A_x .* conj(carriers(:, n)); etc.
 129   r_x = rxFilter(A_x, rolloff, span, sps);
 130   r_y = rxFilter(A_y, rolloff, span, sps);
 131   % Rx filter performs downsampling as well, keep track of this
 132   sps = 2;
 133   Tsamp = Tsamp * spsOrig / sps;
 134
 135   %% Rx shot noise
 136   photonpersym = mean(abs(r_x) .^ 2) / Rsym / Eperphoton;
 137   snr = photonpersym;
 138   r_x = awgn(r_x, snr, 'measured', 'linear');
 139   r_y = awgn(r_y, snr, 'measured', 'linear');
 140
 141   %% -- Begin DSP channel equalization --
 142   %% Chromatic dispersion compensation
 143   r_x = CDCompensation(r_x, D, lambda, z, Tsamp);
 144   r_y = CDCompensation(r_y, D, lambda, z, Tsamp);
 145   r_x = r_x(2:2:end);
 146   r_y = r_y(2:2:end);
 147
 148   %% Adaptive filter
 149   [r_x, r_y] = pdm_adaptiveCMA(r_x, r_y);
 150
 151   %% Phase noise correction
 152   r_x = phaseNoiseCorr(r_x, M, 0, 40).';
 153   r_y = phaseNoiseCorr(r_y, M, 0, 40).';
 154
 155   %% Demodulate DE-QPSK
 156   demod_x = deqpskdemod(r_x);
 157   demod_y = deqpskdemod(r_y);
 158
 159   %% Calculate and store BER
 160   [~, ber(i)] = biterr([data_x(:, 1); data_y(:, 1)], [demod_x; demod_y]);
 161
 162   q = 20 * log10(erfcinv(2*ber)*sqrt(2));
 163   if plotlen > 1
 164     figure(fig);
 165     plot(zs, q);
 166   end
 167 end
 168
 169 ber
 170 q