function xCDKerr = splitstepfourier(x, D, lambda, z, Tsamp, gamma)
  %% Simulate chromatic dispersion and Kerr effect,
  %% with attenuation and amplification.
  %% Params:
  %%  - x: input waveform (pulse-shaped)
  %%  - D: dispersion coefficient (ps / (nm km))
  %%  - lambda: wavelength (nm)
  %%  - z: length of fibre (km)
  %%  - Tsamp: sampling time (s)
  %% Output:
  %%  - xCDKerr: x after being dispersed.

  %% Convert everything to SI base units
  c = 299792458; % m/s
  D = D * 1e-6; % s/m^2
  lambda = lambda * 1e-9; % m
  z = z * 1e3; % m

  gamma = gamma * 1e-3; % watt^-1 / m
  dz = 1; % km
  dz = dz * 1e3; % m
  stepnum = z / dz;

  alpha = 0.2; % dB/km
  alpha = alpha / 10 * log(10); % /km
  alpha = alpha * 1e-3; % /m

  hhz = 1j * gamma * dz;


  xCDKerr = x .* exp(abs(x) .^ 2 .* hhz / 2 - alpha * dz / 4);

  for i = 1 : stepnum
    xDFT = fft(xCDKerr);
    n = length(xCDKerr);
    fs = 1 / Tsamp;

    omega = (2*pi * fs / n * [(0 : floor((n-1)/2)), (-ceil((n-1)/2) : -1)]).';
    dispDFT = xDFT .* exp(-1j * omega.^2 * D * lambda^2 * dz / (4 * pi * c));
    xCDKerr = ifft(dispDFT);

    xCDKerr = xCDKerr .* exp(abs(xCDKerr) .^ 2 .* hhz - alpha * dz / 2);
    if mod(i, 50) == 0
      xCDKerr = xCDKerr * 10; % amplification
    end
  end
  xCDKerr = xCDKerr .* exp(-abs(xCDKerr) .^ 2 .* hhz / 2 + alpha * dz / 4);
end
%% References
%% [1]: S.J. Savory, Digital filters for coherent optical receivers, 2008.