Added technical milestone report and changes to 1st presentation

[4yp.git] / chromaticDispersion_FFT.m

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.