Working Kerr effect; PDM; speedups; removed unused files

[4yp.git] / splitstepfourier.m

diff --git a/splitstepfourier.m b/splitstepfourier.m
index f71828d..f0ec5fb 100644 (file)
--- a/splitstepfourier.m
+++ b/splitstepfourier.m
@@ -1,5 +1,6 @@
-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))
@@ -7,8 +8,7 @@ function [xCD, kstart] = splitstepfourier(x, D, lambda, z, Tsamp,gamma,stepnum)
   %%  - 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
@@ -17,40 +17,34 @@ function [xCD, kstart] = splitstepfourier(x, D, lambda, z, Tsamp,gamma,stepnum)
   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.