Added presentation; DE-QPSK; CD with FFT; split-step Fourier

[4yp.git] / chromaticDispersion_FFT.m

diff --git a/chromaticDispersion_FFT.m b/chromaticDispersion_FFT.m
new file mode 100644 (file)
index 0000000..b091485
--- /dev/null
+++ b/chromaticDispersion_FFT.m
@@ -0,0 +1,43 @@
+function [xCD, kstart] = chromaticDispersion_FFT(x, D, lambda, z, Tsamp)
+  %% Simulate chromatic dispersion.
+  %% 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:
+  %%  - xCD: x after being dispersed. Energy of xCD is not normalized.
+  %%  - kstart: starting index of the discrete signal
+
+  %% 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
+
+  %% 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
+
+  xDFT = fft(x);
+  n = length(x);
+  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 * z / (4 * pi * c));
+
+  xCD = ifft(dispDFT);
+
+  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)];
+  end
+  %% pad zeros for compatibility with time-domain filter
+  xCD = [zeros(floor((kmax-1)/2), 1); xCD; zeros(ceil((kmax+1)/2), 1)];
+end
+%% References
+%% [1]: S.J. Savory, Digital filters for coherent optical receivers, 2008.