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