splitstepfourier.m

   1 function [xCD, kstart] = splitstepfourier(x, D, lambda, z, Tsamp,gamma,stepnum)
   2   %% Simulate chromatic dispersion.
   3   %% Params:
   4   %%  - x: input waveform (pulse-shaped)
   5   %%  - D: dispersion coefficient (ps / (nm km))
   6   %%  - lambda: wavelength (nm)
   7   %%  - z: length of fibre (km)
   8   %%  - Tsamp: sampling time (s)
   9   %% Output:
  10   %%  - xCD: x after being dispersed. Energy of xCD is not normalized.
  11   %%  - kstart: starting index of the discrete signal
  12
  13   %% Convert everything to SI base units
  14   c = 299792458; % m/s
  15   D = D * 1e-6; % s/m^2
  16   lambda = lambda * 1e-9; % m
  17   z = z * 1e3; % m
  18
  19   gamma = gamma * 1e-3; % watt^-1 / m
  20   %stepnum = 1
  21   dz = z / stepnum;
  22
  23   hhz = 1j * gamma * dz;
  24
  25
  26   %% Time domain filter length, needed for compatibility with
  27   %% time-domain (convolution) chromatic dispersion.
  28   kmax = floor(abs(D) * lambda^2 * z / (2 * c * Tsamp^2));
  29   kstart = 1 - kmax;
  30
  31   x = [zeros(kmax, 1); x]; % prepend zeros to allow for transients
  32
  33   xCD = x .* exp(abs(x) .^ 2 .* hhz / 2);
  34   for i = 1 : stepnum
  35     xDFT = fft(xCD);
  36     n = length(xCD);
  37     fs = 1 / Tsamp;
  38
  39     omega = (2*pi * fs / n * [(0 : floor((n-1)/2)), (-ceil((n-1)/2) : -1)]).';
  40     dispDFT = xDFT .* exp(-1j * omega.^2 * D * lambda^2 * dz / (4 * pi * c));
  41
  42     xCD = ifft(dispDFT);
  43
  44     xCD = xCD .* exp(abs(xCD) .^ 2 .* hhz / 2);
  45   end
  46   xCD = xCD .* exp(-abs(xCD) .^ 2 .* hhz / 2);
  47
  48   if ceil(kmax/2) > 0
  49     %% fix the order of the samples due to prepending zeros before FFT
  50     xCD = [xCD(ceil(kmax/2):end); xCD(1:ceil(kmax/2)-1)];
  51   end
  52   %% pad zeros for compatibility with time-domain filter
  53   xCD = [zeros(floor((kmax-1)/2), 1); xCD; zeros(ceil((kmax+1)/2), 1)];
  54 end
  55 %% References
  56 %% [1]: S.J. Savory, Digital filters for coherent optical receivers, 2008.