1 function [xCD, kstart] = splitstepfourier(x, D, lambda, z, Tsamp,gamma,stepnum)
2 %% Simulate chromatic dispersion.
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)
10 %% - xCD: x after being dispersed. Energy of xCD is not normalized.
11 %% - kstart: starting index of the discrete signal
13 %% Convert everything to SI base units
16 lambda = lambda * 1e-9; % m
19 gamma = gamma * 1e-3; % watt^-1 / m
23 hhz = 1j * gamma * dz;
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));
31 x = [zeros(kmax, 1); x]; % prepend zeros to allow for transients
33 xCD = x .* exp(abs(x) .^ 2 .* hhz / 2);
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));
44 xCD = xCD .* exp(abs(xCD) .^ 2 .* hhz / 2);
46 xCD = xCD .* exp(-abs(xCD) .^ 2 .* hhz / 2);
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)];
52 %% pad zeros for compatibility with time-domain filter
53 xCD = [zeros(floor((kmax-1)/2), 1); xCD; zeros(ceil((kmax+1)/2), 1)];
56 %% [1]: S.J. Savory, Digital filters for coherent optical receivers, 2008.