-function [xCD, kstart] = splitstepfourier(x, D, lambda, z, Tsamp,gamma,stepnum)
- %% Simulate chromatic dispersion.
+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))
%% - z: length of fibre (km)
%% - Tsamp: sampling time (s)
%% Output:
- %% - xCD: x after being dispersed. Energy of xCD is not normalized.
- %% - kstart: starting index of the discrete signal
+ %% - xCDKerr: x after being dispersed.
%% Convert everything to SI base units
c = 299792458; % m/s
z = z * 1e3; % m
gamma = gamma * 1e-3; % watt^-1 / m
- %stepnum = 1
- dz = z / stepnum;
+ dz = 1; % km
+ dz = dz * 1e3; % m
+ stepnum = z / dz;
- hhz = 1j * gamma * dz;
+ alpha = 0.2; % dB/km
+ alpha = alpha / 10 * log(10); % /km
+ alpha = alpha * 1e-3; % /m
+ hhz = 1j * gamma * dz;
- %% Time domain filter length, needed for compatibility with
- %% time-domain (convolution) chromatic dispersion.
- kmax = floor(abs(D) * lambda^2 * z / (2 * c * Tsamp^2));
- kstart = 1 - kmax;
- x = [zeros(kmax, 1); x]; % prepend zeros to allow for transients
+ xCDKerr = x .* exp(abs(x) .^ 2 .* hhz / 2 - alpha * dz / 4);
- xCD = x .* exp(abs(x) .^ 2 .* hhz / 2);
for i = 1 : stepnum
- xDFT = fft(xCD);
- n = length(xCD);
+ 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);
- xCD = ifft(dispDFT);
-
- xCD = xCD .* exp(abs(xCD) .^ 2 .* hhz / 2);
- end
- xCD = xCD .* exp(-abs(xCD) .^ 2 .* hhz / 2);
-
- if ceil(kmax/2) > 0
- %% fix the order of the samples due to prepending zeros before FFT
- xCD = [xCD(ceil(kmax/2):end); xCD(1:ceil(kmax/2)-1)];
+ xCDKerr = xCDKerr .* exp(abs(xCDKerr) .^ 2 .* hhz - alpha * dz / 2);
+ if mod(i, 50) == 0
+ xCDKerr = xCDKerr * 10; % amplification
+ end
end
- %% pad zeros for compatibility with time-domain filter
- xCD = [zeros(floor((kmax-1)/2), 1); xCD; zeros(ceil((kmax+1)/2), 1)];
+ xCDKerr = xCDKerr .* exp(-abs(xCDKerr) .^ 2 .* hhz / 2 + alpha * dz / 4);
end
%% References
%% [1]: S.J. Savory, Digital filters for coherent optical receivers, 2008.