From: Adrian Iain Lam <adrianiainlam@users.noreply.github.com>
Date: Tue, 12 Feb 2019 01:37:28 +0000 (+0000)
Subject: Working Kerr effect; PDM; speedups; removed unused files
X-Git-Url: https://adrianiainlam.tk/git/?a=commitdiff_plain;h=HEAD;p=4yp.git

Working Kerr effect; PDM; speedups; removed unused files

Kerr effect is now working, showing the expected curve shape as Tx
power varies. Currently this uses a fixed step size of 1km, with
attenuation, and amplification every 50km. Amplification will be removed
later to simulate a true passive network, but now kept to verify results.

Since non-linearity is now introduced, we change the Tx signal to
having 8 samples/sec (previously 2). The Rx then downsamples at the
Rx filter.

"channel.m" is now added to simulate the overall channel with all
effects considered. Currently this includes Kerr, CD and phase
noise. PDM was then implemented, with "pdmchannel.m" being the
PDM analogue of "channel.m", with corresponding new files for the
PDM version of the CMA equalizer and the split-step solver.

Simulations are now significantly faster by doing all filtering
in the frequency domain (so multiplication instead of convolution).
The behaviour is changed slightly: instead of zero initial conditions,
we now take periodic (circular) boundary conditions. Simulation
lengths were also changed to an integer power of 2 to have the
most efficient FFT.

Some initial testing files that are no longer needed are removed.
---

diff --git a/CDCompensation.m b/CDCompensation.m
index 1f3cc2c..9f01950 100644
--- a/CDCompensation.m
+++ b/CDCompensation.m
@@ -21,11 +21,20 @@ function [yCDComp, kstart] = CDCompensation(y, D, lambda, z, Tsamp)
   kstart = -floor(N/2);
 
   % h: FIR filter
-  h = sqrt(j * c * Tsamp^2 / (D * lambda^2 * z)) * ...
-      exp(-j * pi * c * Tsamp^2 * k .^ 2 / (D * lambda^2 * z));
+  h = exp(-j * pi * c * Tsamp^2 * k .^ 2 / (D * lambda^2 * z));
+
+  len_fft = max(length(y), length(h));
+  H = fft(h, len_fft);
+  Y = fft(y, len_fft);
+
+  yCDComp = ifft(H.' .* Y);
+  l = (length(h) - 1) / 2;
+  if l > 0
+    yCDComp = [yCDComp(l:end); yCDComp(1:l-1)];
+  else
+    yCDComp = [yCDComp(end); yCDComp(1:end-1)];
+  end
 
-  yCDComp = upfirdn(y, h);
-  yCDComp = yCDComp(N:end); % truncate filter transients
 end
 %% References
 %% [1]: S.J. Savory, Digital filters for coherent optical receivers, 2008.
diff --git a/CD_AWGN.m b/CD_AWGN.m
index 4db7aad..241d22a 100644
--- a/CD_AWGN.m
+++ b/CD_AWGN.m
@@ -1,18 +1,18 @@
-numSymbs = 5e5;
+numSymbs = 2^16;
 M = 4;
 
 Rsym = 2.5e10; % symbol rate (sym/sec)
 
 rolloff = 0.25;
 span = 6; % filter span
-sps = 2; % samples per symbol
+sps = 8; % samples per symbol
 
 fs = Rsym * sps; % sampling freq (Hz)
 Tsamp = 1 / fs;
 
 t = (0 : 1 / fs : numSymbs / Rsym + (1.5 * span * sps - 1) / fs).';
 
-EbN0_db = 0:0.2:14;
+EbN0_db = 0:0.5:14;
 EbN0 = 10 .^ (EbN0_db ./ 10);
 
 Es = 1;
@@ -36,26 +36,30 @@ x = txFilter(modData, rolloff, span, sps);
 %% Simulate chromatic dispersion
 D = 17; % ps / (nm km)
 lambda = 1550; % nm
-z = 60;%000; % km
+z = 3000; % km
 
-usingFFT = 1
-xCD = chromaticDispersion_FFT(x, D, lambda, z, Tsamp);
-%%xCD = normalizeEnergy(xCD, numSymbs, 1);
-%%xCD = x;
+
+[xCD, xCDkstart] = chromaticDispersion(x, D, lambda, z, Tsamp);
+
+TsampOrig = Tsamp;
 
 for i = 1:plotlen
+  sps = 8;
+
   snr = EbN0_db(i) + 10 * log10(log2(M)) - 10 * log10(sps);
-  noiseEnergy = 10 ^ (-snr / 10);
 
   y = awgn(xCD, snr, 'measured');
-  %%y = xCD;
 
   r = rxFilter(y, rolloff, span, sps);
-  rCDComp = CDCompensation(r, D, lambda, z, Tsamp);
+
+  sps = 2;
+  Tsamp = TsampOrig * 4;
+
+  [rCDComp, CDCompkstart] = CDCompensation(r, D, lambda, z, Tsamp);
   rCDComp = normalizeEnergy(rCDComp, numSymbs*sps, 1);
 
-  rSampled = rCDComp(sps*span/2+1:sps:(numSymbs+span/2)*sps);
-  rNoCompSampled = r(sps*span/2+1:sps:(numSymbs+span/2)*sps);
+  rSampled = rCDComp(2:2:end);
+  rNoCompSampled = r(2:2:end);
 
   %% rotate rNoCompSampled to match original data
   theta = angle(-sum(rNoCompSampled .^ M)) / M;
diff --git a/adaptiveCMA.m b/adaptiveCMA.m
index 084eab4..41b1e48 100644
--- a/adaptiveCMA.m
+++ b/adaptiveCMA.m
@@ -1,4 +1,4 @@
-function adaptFilterOut = adaptiveCMA(rSampled)
+function [adaptFilterOut, convergeIdx] = adaptiveCMA(rSampled)
   %% adaptive filter
   %% CMA
   taps = 19; % ODD taps
@@ -17,6 +17,7 @@ function adaptFilterOut = adaptiveCMA(rSampled)
   adaptFilterOut = zeros(numSymbs, 1);
   converged = 0;
   convergeCount = 0;
+  convergeIdx = Inf;
 
   for it = 1:numSymbs
     if it <= (taps - 1) / 2;
@@ -30,14 +31,14 @@ function adaptFilterOut = adaptiveCMA(rSampled)
     xout = sum(hxx .* xp);
     ex = 1 - abs(xout) ^ 2;
 
-    if abs(ex) < 1e-3
+    if abs(ex) < 0.01
       convergeCount = convergeCount + 1;
     else
       convergeCount = 0;
     end
     if ~converged && convergeCount >= 10
-      converged = 1
-      it
+      converged = 1;
+      convergeIdx = it;
     end
 
     adaptFilterOut(it) = xout;
diff --git a/agrawalAppendixB.m b/agrawalAppendixB.m
deleted file mode 100644
index 295e508..0000000
--- a/agrawalAppendixB.m
+++ /dev/null
@@ -1,58 +0,0 @@
-distance = input('Enter fiber length in L_D ');
-beta2 = input('dispersion: 1 for normal, -1 for anomalous ');
-N = input('Nonlinear parameter N = '); % Soliton order
-mshape = input('m = 0 for sech, m > 0 for super-Gaussian ');
-chirp0 = 0;
-
-% set simulation parameters
-nt = 1024; Tmax = 32; % FFT points and window size
-step_num = round(20 * distance * N^2); % No. of z steps
-deltaz = distance/step_num; % step size in z
-dtau = (2*Tmax) / nt; % step size in tau
-
-%% tau and omega arrays
-tau = (-nt/2 : nt/2-1) * dtau; % temporal grid
-omega = (pi/Tmax) * [(0:nt/2-1) (-nt/2:-1)]; % freq grid
-
-if mshape == 0
-  uu = sech(tau) .* exp(-0.5j * chirp0 * tau.^2);
-else
-  uu = exp(-0.5 * (1 + 1j * chirp0) .* tau.^(2 * mshape));
-end
-
-%% plot input pulse shape and spectrum
-temp = fftshift(ifft(uu)) .* (nt * dtau) / sqrt(2 * pi); % spectrum
-figure(1); clf; subplot(2,1,1);
-plot(tau, abs(uu).^2, '--k'); hold on;
-axis([-5 5 0 inf]);
-xlabel('Normalized Time');
-ylabel('Normalized Power');
-title('Input and Output pulse shape and spectrum');
-
-subplot(2, 1, 2);
-plot(fftshift(omega)/(2*pi), abs(temp) .^ 2, '--k'); hold on;
-axis([-.5 .5 0 inf]);
-xlabel('Normaized freq');
-ylabel('spectral power');
-
-%% store dispersive phase shifts to speed up code
-dispersion = exp(0.5j * beta2 * omega.^2 * deltaz); % [hase factor
-hhz = 1j * N^2 * deltaz;
-
-%% begin main loop
-%% N/2 -> D -> N/2 first half step nonlinear
-temp = uu .* exp(abs(uu) .^ 2 .* hhz / 2);
-for n = 1 : step_num
-  f_temp = ifft(temp) .* dispersion;
-  uu = fft(f_temp);
-  temp = uu .* exp(abs(uu) .^ 2 .* hhz);
-end
-uu = temp .* exp(-abs(uu) .^ 2 .* hhz); % final field
-temp = fftshift(ifft(uu)) .* (nt*dtau) / sqrt(2 * pi); % final spectrum
-%% end of main loop
-
-%% plot output pulse shape and spectrum
-subplot(2,1,1);
-plot(tau, abs(uu).^2, '-k');
-subplot(2,1,2);
-plot(fftshift(omega) / (2*pi), abs(temp).^2, '-k');
diff --git a/baseband.m b/baseband.m
deleted file mode 100644
index d6f0f1b..0000000
--- a/baseband.m
+++ /dev/null
@@ -1,69 +0,0 @@
-function baseband(rolloff, M, numSymbs)
-  %% Set defaults for inputs
-  if nargin < 3
-    numSymbs = 1000;
-  end
-  if nargin < 2
-    M = 2;
-  end
-  if nargin < 1
-    rolloff = 0.5;
-  end
-
-  %% https://www.mathworks.com/help/comm/examples/passband-modulation-with-adjacent-channel-interference.html
-  Rsym = 1e6; % symbol rate (sym/sec)
-
-  span = 6; % filter span
-  sps = 2; % samples per symbol
-
-  fs = Rsym * sps; % sampling freq (Hz)
-
-  t = (0 : 1 / fs : numSymbs / Rsym + (1.5 * span * sps - 1) / fs)';
-
-  EbN0_db = 0:0.2:10;
-  EbN0 = 10 .^ (EbN0_db ./ 10);
-  Es = 1;
-  Eb = Es / log2(M);
-  N0 = Eb ./ EbN0;
-
-  EsN0 = EbN0 .* log2(M);
-  EsN0_db = 10 .* log10(EsN0);
-
-  plotlen = length(EbN0);
-  ber = zeros(1, plotlen);
-
-  data = randi([0 M - 1], numSymbs, 1);
-  modData = pskmod(data, M, 0, 'gray');
-
-  xBaseband = txFilter(modData, rolloff, span, sps);
-
-  for i = 1:plotlen
-    snr = EbN0_db(i) + 10 * log10(log2(M)) - 10 * log10(sps);
-    noiseEnergy = 10 ^ (-snr / 10);
-
-    yBaseband = awgn(xBaseband, snr, 'measured');
-
-    rBaseband = rxFilter(yBaseband, rolloff, span, sps);
-
-    rSampled = rBaseband(sps*span/2+1:sps:(numSymbs+span/2)*sps);
-    demodData = pskdemod(rSampled, M, 0, 'gray');
-    [bitErrors, ber(i)] = biterr(data, demodData);
-  end
-
-  figure(1);
-  clf;
-
-  %% Plot simulated results
-  semilogy(EbN0_db, ber, 'r', 'LineWidth', 2);
-  hold on;
-
-  theoreticalPSK(EbN0_db, M, 'b', 'LineWidth', 1);
-  legend('Simulated', 'Discrete');
-
-  title(strcat(num2str(M), '-PSK with Gray code'));
-  grid on;
-  xlabel('$E_b/N_0$ (dB)');
-  ylabel('BER');
-
-  formatFigure;
-end
diff --git a/channel.m b/channel.m
new file mode 100644
index 0000000..260d21f
--- /dev/null
+++ b/channel.m
@@ -0,0 +1,123 @@
+numSymbs = 2^16;
+M = 4;
+
+Rsym = 2.5e10; % symbol rate (sym/sec)
+Tsym = 1 / Rsym; % symbol period (sec)
+
+rolloff = 0.25;
+span = 6; % filter span
+sps = 8; % samples per symbol
+
+fs = Rsym * sps; % sampling freq (Hz)
+Tsamp = 1 / fs;
+
+t = (0 : 1 / fs : numSymbs / Rsym + (1.5 * span * sps - 1) / fs).';
+
+power_dBm = -6:1:4;
+%%power_dBm = 0;
+power = 10 .^ (power_dBm / 10) * 1e-3; % watts
+
+Es = power * Tsym; % joules
+Eb = Es / log2(M); % joules
+
+N0ref_db = 10; % Eb/N0 at power = 1mW
+%% Fix N0, such that Eb/N0 = N0ref_db at power = 1mW
+N0 = 1e-3 * Tsym / (log2(M) * 10 ^ (N0ref_db / 10)); % joules
+%% At current settings, N0 = 0.002 pJ
+
+plotlen = length(power);
+
+ber = zeros(1, plotlen);
+
+data = randi([0 M - 1], numSymbs, 1);
+%%modData = dpskmod(data, M, 0, 'gray');
+modData = pskmod(data, M, 0, 'gray');
+for i = 2:numSymbs
+  modData(i) = modData(i) * modData(i-1);
+end
+
+
+%% Chromatic dispersion
+D = 17; % ps / (nm km)
+lambda = 1550; % nm
+z = 100; % km
+
+
+linewidthTx = 0; % Hz
+linewidthLO = 1e6; % Hz
+
+
+TsampOrig = Tsamp;
+
+x_P1 = txFilter(modData, rolloff, span, sps);
+
+
+for i = 1:plotlen
+  sps = 8;
+  Tsamp = TsampOrig;
+
+  snr = Es(i) / sps / N0;
+  snr_dB = 10 * log10(snr);
+
+  %%x = txFilter(modData, rolloff, span, sps);
+  %% Now, sum(abs(x) .^ 2) / length(x) should be 1.
+  %% We can set its power simply by multiplying.
+  x = sqrt(power(i)) * x_P1;
+
+  %% We can now do split-step Fourier.
+  gamma = 1.2; % watt^-1 / km
+
+
+  xCDKerr = splitstepfourier(x, D, lambda, z, Tsamp, gamma);
+
+  xpn = phaseNoise(xCDKerr, linewidthTx, linewidthLO, Tsamp);
+
+  y = awgn(xpn, snr_dB, 'measured', 'db');
+  %y = xCDKerr;
+
+  r = rxFilter(y, rolloff, span, sps);
+  sps = 2;
+  Tsamp = Tsamp * 4;
+
+  rCDComp = CDCompensation(r, D, lambda, z, Tsamp);
+  rCDComp = normalizeEnergy(rCDComp, numSymbs * sps, 1);
+
+  rSampled = rCDComp(2:2:end);
+
+  %% adaptive filter
+  [adaptFilterOut, convergeIdx] = adaptiveCMA(rSampled);
+
+  pncorr = phaseNoiseCorr(adaptFilterOut, M, 0, 40).';
+
+  demodAdapt = pskdemod(pncorr, M, 0, 'gray');
+  remod = pskmod(demodAdapt, M, 0, 'gray');
+  delayed = [1; remod(1:end-1)];
+  demod = pskdemod(remod .* conj(delayed), M, 0, 'gray');
+
+  if convergeIdx < Inf
+    [~, ber(i)] = biterr(data(convergeIdx:end), demod(convergeIdx:end));
+  else
+    [~, ber(i)] = biterr...
+                    (data(ceil(0.8*numSymbs):end), ...
+                     demod(ceil(0.8*numSymbs):end));
+  end
+end
+
+ber
+
+
+figure(1);
+clf;
+
+%% Plot simulated results
+qp = 20 * log10(erfcinv(2*ber)*sqrt(2));
+plot(power_dBm, qp, 'Color', [0, 0.6, 0], 'LineWidth', 2);
+hold on;
+
+title({'CD + Kerr + CD compensation', ...
+       strcat(['$D = 17$ ps/(nm km), $z = ', num2str(z), '$ km'])});
+grid on;
+xlabel('Optical power (dBm)');
+ylabel('$20 \log_{10}\left(\sqrt{2}\mathrm{erfc}^{-1}(2 BER)\right)$');
+
+formatFigure;
diff --git a/chromaticDispersion.m b/chromaticDispersion.m
index 49a49a3..bc04e59 100644
--- a/chromaticDispersion.m
+++ b/chromaticDispersion.m
@@ -26,10 +26,21 @@ function [xCD, kstart] = chromaticDispersion(x, D, lambda, z, Tsamp)
   t = k * Tsamp;
 
   % Impulse response
-  g = sqrt(c / (j * D * lambda^2 * z)) * ...
-      exp(j * pi * c / (D * lambda^2 * z) * t .^ 2);
+  g = exp(j * pi * c / (D * lambda^2 * z) * t .^ 2);
 
-  xCD = conv(x, g);
+  lenx = length(x);
+  leng = length(g);
+
+  len_fft = max(lenx, leng);
+
+  G = fft(g, len_fft);
+  X = fft(x, len_fft);
+
+  xCD = ifft(G.' .* X);
+  l = (leng - 1) / 2;
+  if l > 0
+    xCD = [xCD(l:end); xCD(1:l-1)];
+  end
 
   kstart = 1 - kmax;
 end
diff --git a/chromaticDispersion1Signal.m b/chromaticDispersion1Signal.m
index dd169ac..6441dd5 100644
--- a/chromaticDispersion1Signal.m
+++ b/chromaticDispersion1Signal.m
@@ -5,7 +5,7 @@ Rsym = 2.5e10; % symbol rate (sym/sec)
 
 span = 6; % Tx/Rx filter span
 rolloff = 0.25; % Tx/Rx RRC rolloff
-sps = 2; % samples per symbol
+sps = 8; % samples per symbol
 
 fs = Rsym * sps; % sampling freq (Hz)
 Tsamp = 1 / fs;
@@ -21,9 +21,9 @@ x = normalizeEnergy(x, numSymbs*sps, 1);
 %% Simulate chromatic dispersion
 D = 17; % ps / (nm km)
 lambda = 1550; % nm
-z = 5000 % km
+z = 500 % km
 
-[xCD, xCDkstart] = chromaticDispersion_FFT(x, D, lambda, z, Tsamp);
+[xCD, xCDkstart] = chromaticDispersion(x, D, lambda, z, Tsamp);
 
 EbN0_db = 8;
 snr = EbN0_db + 10 * log10(log2(M)) - 10 * log10(sps);
@@ -32,11 +32,16 @@ snr = EbN0_db + 10 * log10(log2(M)) - 10 * log10(sps);
 y = xCD;
 
 r = rxFilter(y, rolloff, span, sps);
+
+sps = 2;
+Tsamp = Tsamp * 4;
+
+
 [rCDComp, CDCompkstart] = CDCompensation(r, D, lambda, z, Tsamp);
 rCDComp = normalizeEnergy(rCDComp, numSymbs*sps, 1);
 
-rSampled = rCDComp(sps*span/2+1:sps:(numSymbs+span/2)*sps);
-rNoCompSa = r(sps*span/2+1:sps:(numSymbs+span/2)*sps);
+rSampled = rCDComp(2:2:end);
+rNoCompSa = r(2:2:end);
 
 %% if no CD comp, then rotate constellation. Use:
 theta = angle(-sum(rNoCompSa .^ M)) / M;
@@ -57,7 +62,6 @@ end
 rSampled = rSampled .* exp(-1j * theta);
 
 
-
 %%rAdaptEq = adaptiveCMA(rSampled);
 %{
 %% Compare original signal and compensated signal
diff --git a/chromaticDispersionTest.m b/chromaticDispersionTest.m
deleted file mode 100644
index 513738f..0000000
--- a/chromaticDispersionTest.m
+++ /dev/null
@@ -1,35 +0,0 @@
-Rsym = 2.5e10; % symbol/sec
-sps = 4; % samples/symbol
-fs = Rsym * sps; % sampling freq (Hz)
-Tsamp = 1 / fs;
-
-x = 10 * normpdf(1:30, 10, 3);
-
-xCD = chromaticDispersion(x, 20, 1550, 30, Tsamp);
-
-figure(10001);
-subplot(121);
-plot(abs(x) .^ 2);
-axis([0 30 -Inf Inf]);
-title('Before dispersion');
-ylabel('Power (arb. unit)');
-xlabel('Time (arb. unit)');
-subplot(122);
-plot(abs(xCD) .^ 2)
-axis([20 50 -Inf Inf]);
-title('After dispersion');
-xlabel('Time (arb.~unit)');
-ylabel('Power (arb.~unit)');
-
-subplot(121);
-formatFigure;
-subplot(122)
-formatFigure;
-set(figure(10001), 'Units', 'centimeters', ...
-    'Position', [0 0 24 8], 'PaperPositionMode', 'auto');
-subplot(121);
-set(gca, 'XTick', 0:5:30);
-subplot(122);
-set(gca, 'XTick', 20:5:50);
-
-%saveas(gcf, 'dispersionTest.png');
diff --git a/chromaticDispersion_FFT.m b/chromaticDispersion_FFT.m
deleted file mode 100644
index b091485..0000000
--- a/chromaticDispersion_FFT.m
+++ /dev/null
@@ -1,43 +0,0 @@
-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.
diff --git a/discretePSK_BER_SNR.m b/discretePSK_BER_SNR.m
deleted file mode 100644
index cdff01f..0000000
--- a/discretePSK_BER_SNR.m
+++ /dev/null
@@ -1,57 +0,0 @@
-function discretePSK_BER_SNR(M, numSymbs)
-  %% Set defaults for inputs
-  if nargin < 2
-    numSymbs = 1000;
-  end
-  if nargin < 1
-    M = 2;
-  end
-
-  EbN0_db = 0:0.2:10;
-  EbN0 = 10 .^ (EbN0_db ./ 10);
-
-  Es = 1;
-  Eb = Es / log2(M);
-  N0 = Eb ./ EbN0;
-
-  EsN0 = EbN0 .* log2(M);
-  EsN0_db = 10 .* log10(EsN0);
-
-  plotlen = length(EbN0);
-
-  ber = zeros(1, plotlen);
-
-  data = randi([0, M - 1], 1, numSymbs);
-  txsig = pskmod(data, M, 0, 'gray');
-
-  for i = 1:plotlen
-    rxsig = awgn(txsig, EsN0_db(i));
-    demodData = pskdemod(rxsig, M, 0, 'gray');
-
-    [bitErrors, ber(i)] = biterr(data, demodData);
-  end
-
-  figure(1);
-  clf;
-
-  %% Plot simulated results
-  semilogy(EbN0_db, ber, 'r', 'LineWidth', 2);
-  hold on;
-
-  theoreticalPSK(EbN0_db, M, 'b', 'LineWidth', 1);
-
-  if M == 2 || M == 4
-    legend('Simulated', 'Theoretical');
-  else
-    legend('Simulated', 'Approximation');
-  end
-
-  title(strcat(num2str(M), '-PSK with Gray code'));
-  grid on;
-  xlabel('$E_b/N_0$ (dB)');
-  ylabel('BER');
-
-  formatFigure;
-  %saveas(gcf, strcat('BER_SNR_', num2str(M), 'PSK_', num2str(numSymbs), ...
-  %                   '.svg'));
-end
diff --git a/formatFigure.m b/formatFigure.m
index a09586d..c3d607f 100644
--- a/formatFigure.m
+++ b/formatFigure.m
@@ -4,3 +4,10 @@ set(gca, 'TickLabelInterpreter', 'Latex');
 set(gca, 'FontSize', 14);
 leg = findobj(gcf, 'Type', 'Legend');
 set(leg, 'FontSize', 14);
+set(leg, 'Interpreter', 'Latex');
+tit = get(gca, 'title');
+set(tit, 'Interpreter', 'Latex');
+xl = get(gca, 'XLabel');
+set(xl, 'Interpreter', 'Latex');
+yl = get(gca, 'YLabel');
+set(yl, 'Interpreter', 'Latex');
diff --git a/kerr.m b/kerr.m
index c67bbb2..7e534dd 100644
--- a/kerr.m
+++ b/kerr.m
@@ -1,4 +1,4 @@
-numSymbs = 5e4;
+numSymbs = 2^16;
 M = 4;
 
 Rsym = 2.5e10; % symbol rate (sym/sec)
@@ -6,40 +6,46 @@ Tsym = 1 / Rsym; % symbol period (sec)
 
 rolloff = 0.25;
 span = 6; % filter span
-sps = 2; % samples per symbol
+sps = 8; % samples per symbol
 
 fs = Rsym * sps; % sampling freq (Hz)
 Tsamp = 1 / fs;
 
 t = (0 : 1 / fs : numSymbs / Rsym + (1.5 * span * sps - 1) / fs).';
 
-
-power_dBm = -3:0.25:4;
+power_dBm = -6:1:4;
+%%power_dBm = 0;
 power = 10 .^ (power_dBm / 10) * 1e-3; % watts
 
 Es = power * Tsym; % joules
 Eb = Es / log2(M); % joules
 
-N0ref_dB = 10; % Eb/N0 at power = 1mW
-%% Fix N0, such that Eb/N0 = N0ref_dB at power = 1mW
-N0 = 1e-3 * Tsym / (log2(M) * 10 ^ (N0ref_dB / 10)); % joules
-
+N0ref_db = 10; % Eb/N0 at power = 1mW
+%% Fix N0, such that Eb/N0 = N0ref_db at power = 1mW
+N0 = 1e-3 * Tsym / (log2(M) * 10 ^ (N0ref_db / 10)); % joules
+%% At current settings, N0 = 0.002 pJ
 
 plotlen = length(power);
 
 ber = zeros(1, plotlen);
 
 data = randi([0 M - 1], numSymbs, 1);
-modData = pskmod(data, M, pi / M, 'gray');
+modData = dpskmod(data, M, 0, 'gray');
+%%modData = pskmod(data, M, pi/4, 'gray');
 
 
 %% Chromatic dispersion
 D = 17; % ps / (nm km)
 lambda = 1550; % nm
-z = 600; % km
+z = 100; % km
 
 
+TsampOrig = Tsamp;
+
 for i = 1:plotlen
+  sps = 8;
+  Tsamp = TsampOrig;
+
   snr = Es(i) / sps / N0;
   snr_dB = 10 * log10(snr);
 
@@ -50,58 +56,52 @@ for i = 1:plotlen
 
   %% We can now do split-step Fourier.
   gamma = 1.2; % watt^-1 / km
-  %%stepnum = round(40 * z * gamma); % Nonlinear Fiber optics, App B
-  stepnum = 100;
-  xCD = splitstepfourier(x, D, lambda, z, Tsamp, gamma, stepnum);
 
-  y = awgn(xCD, snr, power(i), 'linear');
-  %%y = xCD;
+
+  xCDKerr = splitstepfourier(x, D, lambda, z, Tsamp, gamma);
+
+  y = awgn(xCDKerr, snr_dB, 'measured', 'db');
+  %y = xCDKerr;
 
   r = rxFilter(y, rolloff, span, sps);
+  sps = 2;
+  Tsamp = Tsamp * 4;
+
   rCDComp = CDCompensation(r, D, lambda, z, Tsamp);
-  rCDComp = normalizeEnergy(rCDComp, numSymbs*sps, 1);
+  rCDComp = normalizeEnergy(rCDComp, numSymbs * sps, 1);
 
-  rSampled = rCDComp(sps*span/2+1:sps:(numSymbs+span/2)*sps);
-  rNoCompSampled = r(sps*span/2+1:sps:(numSymbs+span/2)*sps);
+  rSampled = rCDComp(2:2:end);
 
-  %% rotate rNoCompSampled to match original data
-  theta = angle(-sum(rNoCompSampled .^ M)) / M;
-  %% if theta approx +pi/M, wrap to -pi/M
-  if abs(theta - pi / M) / (pi / M) < 0.1
-    theta = -pi / M;
-  end
-  rNoCompSampled = rNoCompSampled .* exp(-j * theta);
+  %% adaptive filter
+  [adaptFilterOut, convergeIdx] = adaptiveCMA(rSampled);
 
+  demod = dpskdemod(adaptFilterOut, M, 0, 'gray');
+  %%demod = pskdemod(adaptFilterOut, M, pi/4, 'gray');
 
-  %% Not entirely sure why, but after using FFT instead of time-domain
-  %% convolution for simulating CD, we now need to do the same rotation
-  %% for rSampled as well, but this time with a positive rotation.
-  theta = angle(-sum(rSampled .^ M)) / M;
-  if abs(theta + pi / M) / (pi / M) < 0.1
-    theta = +pi / M;
+  if convergeIdx < Inf
+    [~, ber(i)] = biterr(data(convergeIdx:end), demod(convergeIdx:end));
+  else
+    [~, ber(i)] = biterr...
+                    (data(ceil(0.8*numSymbs):end), ...
+                     demod(ceil(0.8*numSymbs):end));
   end
-  rSampled = rSampled .* exp(-1j * theta);
-
+end
 
-  %% adaptive filter
-  adaptFilterOut = adaptiveCMA(rSampled);
+ber
 
-  demodAdapt = pskdemod(adaptFilterOut, M, pi / M, 'gray');
-  [~, ber(i)] = biterr(data, demodAdapt);
-end
 
-figure(1);
+figure;
 clf;
 
 %% Plot simulated results
-semilogy(power_dBm, ber, 'Color', [0, 0.6, 0], 'LineWidth', 2);
+qp = 20 * log10(erfcinv(2*ber)*sqrt(2));
+plot(power_dBm, qp, 'Color', [0, 0.6, 0], 'LineWidth', 2);
 hold on;
 
 title({'CD + Kerr + CD compensation', ...
-       strcat(['$D = 17$ ps/(nm km), $z = ', num2str(z), '$ km']), ...
-       strcat(['$E_b/N_0 = ', num2str(N0ref_dB), '$ dB at 1 mW'])});
+       strcat(['$D = 17$ ps/(nm km), $z = ', num2str(z), '$ km'])});
 grid on;
 xlabel('Optical power (dBm)');
-ylabel('BER');
+ylabel('$20 \log_{10}\left(\sqrt{2}\mathrm{erfc}^{-1}(2 BER)\right)$');
 
 formatFigure;
diff --git a/kerr1Signal.m b/kerr1Signal.m
deleted file mode 100644
index f4ea292..0000000
--- a/kerr1Signal.m
+++ /dev/null
@@ -1,110 +0,0 @@
-numSymbs = 5e5;
-M = 4;
-
-Rsym = 2.5e10; % symbol rate (sym/sec)
-Tsym = 1 / Rsym; % symbol period (sec)
-
-rolloff = 0.25;
-span = 6; % filter span
-sps = 2; % samples per symbol
-
-fs = Rsym * sps; % sampling freq (Hz)
-Tsamp = 1 / fs;
-
-t = (0 : 1 / fs : numSymbs / Rsym + (1.5 * span * sps - 1) / fs).';
-
-
-%%power_dBm = -3:0.2:4;
-power_dBm = [0];
-power = 10 .^ (power_dBm / 10) * 1e-3; % watts
-
-Es = power * Tsym; % joules
-Eb = Es / log2(M); % joules
-
-N0ref_db = 10; % Eb/N0 at power = 1mW
-%% Fix N0, such that Eb/N0 = N0ref_db at power = 1mW
-N0 = 1e-3 * Tsym / (log2(M) * 10 ^ (N0ref_db / 10)); % joules
-
-
-plotlen = length(power);
-
-ber = zeros(1, plotlen);
-
-data = randi([0 M - 1], numSymbs, 1);
-modData = pskmod(data, M, pi / M, 'gray');
-
-
-%% Chromatic dispersion
-D = 17; % ps / (nm km)
-lambda = 1550; % nm
-z = 600; % km
-
-
-for i = 1:plotlen
-  snr = Es(i) / sps / N0;
-  snr_dB = 10 * log10(snr);
-
-  x = txFilter(modData, rolloff, span, sps);
-  %% Now, sum(abs(x) .^ 2) / length(x) should be 1.
-  %% We can set its power simply by multiplying.
-  x = sqrt(power(i)) * x;
-
-  %% We can now do split-step Fourier.
-  gamma = 1.2; % watt^-1 / km
-  %%stepnum = round(40 * z * gamma); % Nonlinear Fiber optics, App B
-  stepnum = 100;
-  xCD = splitstepfourier(x, D, lambda, z, Tsamp, gamma, stepnum);
-
-  y = awgn(xCD, snr, power(i), 'linear');
-  %%y = xCD;
-
-  r = rxFilter(y, rolloff, span, sps);
-  rCDComp = CDCompensation(r, D, lambda, z, Tsamp);
-  rCDComp = normalizeEnergy(rCDComp, numSymbs*sps, 1);
-
-  rSampled = rCDComp(sps*span/2+1:sps:(numSymbs+span/2)*sps);
-  rNoCompSampled = r(sps*span/2+1:sps:(numSymbs+span/2)*sps);
-
-  %% rotate rNoCompSampled to match original data
-  theta = angle(-sum(rNoCompSampled .^ M)) / M;
-  %% if theta approx +pi/M, wrap to -pi/M
-  if abs(theta - pi / M) / (pi / M) < 0.1
-    theta = -pi / M;
-  end
-  rNoCompSampled = rNoCompSampled .* exp(-j * theta);
-
-
-  %% Not entirely sure why, but after using FFT instead of time-domain
-  %% convolution for simulating CD, we now need to do the same rotation
-  %% for rSampled as well, but this time with a positive rotation.
-  theta = angle(-sum(rSampled .^ M)) / M;
-  if abs(theta + pi / M) / (pi / M) < 0.1
-    theta = +pi / M;
-  end
-  rSampled = rSampled .* exp(-1j * theta);
-
-  %% adaptive filter
-  adaptFilterOut = adaptiveCMA(rSampled);
-
-  demodAdapt = pskdemod(adaptFilterOut, M, pi / M, 'gray');
-  [~, ber(i)] = biterr(data, demodAdapt)
-end
-
-return
-
-
-figure(1);
-clf;
-
-%% Plot simulated results
-semilogy(power_dBm, ber, 'Color', [0, 0.6, 0], 'LineWidth', 2);
-hold on;
-
-title({'CD + Kerr + CD compensation', ...
-       strcat(['$D = 17$ ps/(nm km), $z = ', num2str(z), '$ km'])});
-grid on;
-%%xlabel('$E_b/N_0$ (dB)');
-xlabel('Optical power (dBm)');
-ylabel('BER');
-
-formatFigure;
diff --git a/passband.m b/passband.m
deleted file mode 100644
index 9e88af6..0000000
--- a/passband.m
+++ /dev/null
@@ -1,79 +0,0 @@
-function passband(rolloff, M, numSymbs)
-  %% Set defaults for inputs
-  if nargin < 3
-    numSymbs = 1000;
-  end
-  if nargin < 2
-    M = 2;
-  end
-  if nargin < 1
-    rolloff = 0.5;
-  end
-
-  %% https://www.mathworks.com/help/comm/examples/passband-modulation-with-adjacent-channel-interference.html
-  Rsym = 2.5e10; % symbol rate (sym/sec)
-
-  span = 6; % filter span
-  sps = 4; % samples per symbol
-
-  fs = Rsym * sps; % sampling freq (Hz)
-
-  t = (0 : 1 / fs : numSymbs / Rsym + (1.5 * span * sps - 1) / fs)';
-
-
-  EbN0_db = 0:0.2:10;
-  EbN0 = 10 .^ (EbN0_db ./ 10);
-
-  Es = 1;
-  Eb = Es / log2(M);
-  N0 = Eb ./ EbN0;
-
-  EsN0 = EbN0 .* log2(M);
-  EsN0_db = 10 .* log10(EsN0);
-
-  plotlen = length(EbN0);
-
-  ber = zeros(1, plotlen);
-
-
-  data = randi([0 M - 1], numSymbs, 1);
-  modData = pskmod(data, M, 0, 'gray');
-
-  xBaseband = txFilter(modData, rolloff, span, sps);
-
-  fc = 1e12; % Carrier freq (Hz)
-  carrier = sqrt(2) * exp(j * 2 * pi * fc * t);
-
-  xPassband = xBaseband .* carrier(1:length(xBaseband));
-
-  for i = 1:plotlen
-    snr = EbN0_db(i) + 10 * log10(log2(M)) - 10 * log10(sps);
-    noiseEnergy = 10 ^ (-snr / 10);
-
-    yPassband = awgn(xPassband, snr, 'measured');
-    yLO = yPassband .* conj(carrier(1:length(yPassband)));
-    rBaseband = rxFilter(yLO, rolloff, span, sps);
-
-    rSampled = rBaseband(sps*span/2+1:sps:(numSymbs + span/2) * sps);
-
-    demodData = pskdemod(rSampled, M, 0, 'gray');
-    [bitErrors, ber(i)] = biterr(data, demodData);
-  end
-
-  figure(1);
-  clf;
-
-  %% Plot simulated results
-  semilogy(EbN0_db, ber, 'r', 'LineWidth', 2);
-  hold on;
-
-  theoreticalPSK(EbN0_db, M, 'b', 'LineWidth', 1);
-  legend('Simulated RRC', 'Discrete');
-
-  title(strcat(num2str(M), '-PSK RRC with Gray code'));
-  grid on;
-  xlabel('$E_b/N_0$ (dB)');
-  ylabel('BER');
-
-  formatFigure;
-end
diff --git a/pdm_adaptiveCMA.m b/pdm_adaptiveCMA.m
new file mode 100644
index 0000000..7e6c9eb
--- /dev/null
+++ b/pdm_adaptiveCMA.m
@@ -0,0 +1,72 @@
+function [x, y] = pdm_adaptiveCMA(rx, ry)
+  taps = 19; % ODD taps
+  hxx = zeros(taps, 1);
+  hxy = zeros(taps, 1);
+  hyx = zeros(taps, 1);
+  hyy = zeros(taps, 1);
+  %% hxx: real indices  -K, ..., 0, ..., K.   K = floor(taps/2)
+  %%    MATLAB indices   1      1+K     taps
+  %% initialize hxx, hxx[0] = 1, hxx[k] = hxx[-k] = 0
+  hxx(ceil(taps/2)) = 1;
+  hxy(ceil(taps/2)) = 1;
+  hyx(ceil(taps/2)) = 1;
+  hyy(ceil(taps/2)) = 1;
+
+  mu = 1e-3;
+  numSymbs = length(rx);
+
+  %% Check average energy of symbols
+  rx = rx / sqrt(mean(abs(rx) .^ 2));
+  ry = ry / sqrt(mean(abs(ry) .^ 2));
+
+  x = zeros(numSymbs, 1);
+  y = zeros(numSymbs, 1);
+  converged = 0;
+  convergeCount = 0;
+  convergeIdx = Inf;
+
+  for it = 1:numSymbs
+    if it <= (taps - 1) / 2;
+      xp = [zeros((taps - 1) / 2 - it + 1, 1); rx(1:it + (taps - 1) / 2)];
+      yp = [zeros((taps - 1) / 2 - it + 1, 1); ry(1:it + (taps - 1) / 2)];
+    elseif it + (taps - 1) / 2 > numSymbs
+      xp = [rx(it - (taps - 1) / 2 : end); zeros(it + (taps - 1) / 2 - numSymbs, 1)];
+      yp = [ry(it - (taps - 1) / 2 : end); zeros(it + (taps - 1) / 2 - numSymbs, 1)];
+    else
+      xp = rx(it - (taps - 1) / 2 : it + (taps - 1) / 2);
+      yp = ry(it - (taps - 1) / 2 : it + (taps - 1) / 2);
+    end
+
+    xout = sum(hxx .* xp) + sum(hxy .* yp);
+    yout = sum(hyx .* xp) + sum(hyy .* yp);
+    ex = 1 - abs(xout) ^ 2;
+    ey = 1 - abs(yout) ^ 2;
+
+    if abs(ex) < 0.01
+      convergeCount = convergeCount + 1;
+    else
+      convergeCount = 0;
+    end
+    if ~converged && convergeCount >= 10
+      converged = 1;
+      convergeIdx = it;
+    end
+
+    x(it) = xout;
+    y(it) = yout;
+
+    hxx = hxx + mu * ex * xout * conj(xp);
+    hxy = hxy + mu * ex * xout * conj(yp);
+    hyx = hyx + mu * ey * yout * conj(xp);
+    hyy = hyy + mu * ey * yout * conj(yp);
+
+
+    if sum(abs(hxx - hyx)) < 0.01 && sum(abs(hxy - hyy)) < 0.01
+      hxx = 0.5 * (hxx + flipud(conj(hyy)));
+      hyy = conj(flipud(hxx));
+      hxy = 0.5 * (hxy - conj(flipud(hyx)));
+      hyx = -conj(flipud(hxy));
+    end
+
+  end
+end
diff --git a/pdmchannel.m b/pdmchannel.m
new file mode 100644
index 0000000..997d646
--- /dev/null
+++ b/pdmchannel.m
@@ -0,0 +1,137 @@
+numSymbs = 2^14;
+M = 4;
+
+Rsym = 2.5e10; % symbol rate (sym/sec)
+Tsym = 1 / Rsym; % symbol period (sec)
+
+rolloff = 0.25;
+span = 6; % filter span
+sps = 8; % samples per symbol
+
+fs = Rsym * sps; % sampling freq (Hz)
+Tsamp = 1 / fs;
+
+t = (0 : 1 / fs : numSymbs/2 / Rsym - 1/fs).';
+
+power_dBm = -6:1:4;
+%%power_dBm = 0;
+power = 10 .^ (power_dBm / 10) * 1e-3; % watts
+
+Es = power * Tsym; % joules
+Eb = Es / log2(M); % joules
+
+N0ref_db = 10; % Eb/N0 at power = 1mW
+%% Fix N0, such that Eb/N0 = N0ref_db at power = 1mW
+N0 = 1e-3 * Tsym / (log2(M) * 10 ^ (N0ref_db / 10)); % joules
+%% At current settings, N0 = 0.002 pJ
+
+plotlen = length(power);
+
+ber = zeros(1, plotlen);
+
+data = randi([0 M - 1], numSymbs, 1);
+%%modData = dpskmod(data, M, 0, 'gray');
+modData = pskmod(data, M, 0, 'gray');
+for i = 2:numSymbs
+  modData(i) = modData(i) * modData(i-1);
+end
+
+
+%% Chromatic dispersion
+D = 17; % ps / (nm km)
+lambda = 1550; % nm
+z = 100; % km
+
+linewidthTx = 0; % Hz
+linewidthLO = 1e6; % Hz
+
+TsampOrig = Tsamp;
+
+sig_x = txFilter(modData(1:numSymbs/2), rolloff, span, sps);
+sig_y = txFilter(modData(numSymbs/2+1:end), rolloff, span, sps);
+
+for i = 1:plotlen
+  sps = 8;
+  Tsamp = TsampOrig;
+
+  snr = Es(i) / sps / N0;
+  snr_dB = 10 * log10(snr);
+
+  %%x = txFilter(modData, rolloff, span, sps);
+  %% Now, sum(abs(x) .^ 2) / length(x) should be 1.
+  %% We can set its power simply by multiplying.
+  %%x = sqrt(power(i)) * x;
+  txx = sig_x * sqrt(power(i));
+  txy = sig_y * sqrt(power(i));
+
+  rot_omega = 1e3; % rad/s
+  rot_phi = 2;
+  rot_x = txx .* cos(rot_omega * t) + ...
+          txy .* sin(rot_omega * t) * exp(-1j * rot_phi);
+  rot_y = txx .* -sin(rot_omega * t) * exp(1j * rot_phi) + ...
+          txy .* cos(rot_omega * t);
+
+  %% We can now do split-step Fourier.
+  gamma = 1.2; % watt^-1 / km
+
+  [xCDKerr, yCDKerr] = ssf_pdm(rot_x, rot_y, ...
+                               D, lambda, z, Tsamp, gamma);
+
+  xpn = phaseNoise(xCDKerr, linewidthTx, linewidthLO, Tsamp);
+  ypn = phaseNoise(yCDKerr, linewidthTx, linewidthLO, Tsamp);
+
+  xout = awgn(xpn, snr_dB, 'measured', 'db');
+  yout = awgn(ypn, snr_dB, 'measured', 'db');
+
+  rx = rxFilter(xout, rolloff, span, sps);
+  ry = rxFilter(yout, rolloff, span, sps);
+  sps = 2;
+  Tsamp = Tsamp * 4;
+
+  rxCDComp = CDCompensation(rx, D, lambda, z, Tsamp);
+  ryCDComp = CDCompensation(ry, D, lambda, z, Tsamp);
+
+  rxSampled = rxCDComp(1:2:end);
+  rySampled = ryCDComp(1:2:end);
+
+  %% adaptive filter
+  [xCMA, yCMA] = pdm_adaptiveCMA(rxSampled, rySampled);
+
+  xpncorr = phaseNoiseCorr(xCMA, M, 0, 40).';
+  ypncorr = phaseNoiseCorr(yCMA, M, 0, 40).';
+
+  demodx = pskdemod(xpncorr, M, 0, 'gray');
+  remodx = pskmod(demodx, M, 0, 'gray');
+  delayx = [1; remodx(1:end-1)];
+  demodx = pskdemod(remodx .* conj(delayx), M, 0, 'gray');
+  clear remodx
+  clear delayx
+
+  demody = pskdemod(ypncorr, M, 0, 'gray');
+  remody = pskmod(demody, M, 0, 'gray');
+  delayy = [1; remody(1:end-1)];
+  demody = pskdemod(remody .* conj(delayy), M, 0, 'gray');
+  clear remody
+  clear delayy
+
+
+  [~, ber(i)] = biterr(data, [demodx; demody]);
+end
+
+ber
+
+figure;
+clf;
+
+%% Plot simulated results
+qp = 20 * log10(erfcinv(2*ber)*sqrt(2));
+plot(power_dBm, qp, 'Color', [0, 0.6, 0], 'LineWidth', 2);
+hold on;
+
+title({'CD + Kerr + CD compensation', ...
+       strcat(['$D = 17$ ps/(nm km), $z = ', num2str(z), '$ km'])});
+grid on;
+xlabel('Optical power (dBm)');
+ylabel('$20 \log_{10}\left(\sqrt{2}\mathrm{erfc}^{-1}(2 BER)\right)$');
+
+formatFigure;
diff --git a/phasenoise1signal.m b/phasenoise1signal.m
index f819e65..e2240b9 100644
--- a/phasenoise1signal.m
+++ b/phasenoise1signal.m
@@ -4,7 +4,7 @@ M = 4;
 Rsym = 2.5e10; % symbol rate (sym/sec)
 rolloff = 0.5;
 span = 6; % filter span
-sps = 4; % samples per symbol
+sps = 8; % samples per symbol
 
 fs = Rsym * sps; % sampling freq (Hz)
 Tsamp = 1 / fs;
@@ -39,18 +39,19 @@ linewidthLO = 10e6; % Hz
 
 snr = EbN0_db + 10 * log10(log2(M)) - 10 * log10(sps);
 
-%%y = awgn(xPN, snr, 'measured');
 yPSK = awgn(xPSKpn, snr, 'measured');
 yDPSK = xDPSKpn;
 
 rPSK = rxFilter(yPSK, rolloff, span, sps);
 rDPSK = rxFilter(yDPSK, rolloff, span, sps);
 
+sps = 2;
+Tsamp = Tsamp * 4;
 
-rPSKSa = rPSK(sps*span/2+1:sps:(numSymbs+span/2)*sps);
-rDPSKSa = rDPSK(sps*span/2+1:sps:(numSymbs+span/2)*sps);
+rPSKSa = rPSK(1:2:end);
+rDPSKSa = rDPSK(1:2:end);
 
-[rPSKSaPhEq, phiestsPSK] = phaseNoiseCorr(rPSKSa, M, pi/M, 40);
+[rPSKSaPhEq, phiestsPSK] = phaseNoiseCorr(rPSKSa, M, pi/M, 80);
 
 demodPSK = pskdemod(rPSKSaPhEq, M, pi/M, 'gray');
 demodDPSK = dpskdemod(rDPSKSa, M, pi/M, 'gray');
@@ -60,9 +61,9 @@ demodDPSK = dpskdemod(rDPSKSa, M, pi/M, 'gray');
 
 
 figure(2);
-plot(t(1:40000), repelem(-phiestsPSK, sps));
+plot(t(1:80000), repelem(-phiestsPSK, 8));
 hold on;
-plot(t(1:40000), pTxLoPSK(1:40000));
+plot(t(1:80000), pTxLoPSK(1:80000));
 legend('estimate', 'actual');
 title('Phase noise estimation');
 hold off;
diff --git a/phasenoise_AWGN.m b/phasenoise_AWGN.m
index f055bcc..621a33d 100644
--- a/phasenoise_AWGN.m
+++ b/phasenoise_AWGN.m
@@ -1,11 +1,11 @@
-numSymbs = 5e3;
+numSymbs = 5e4;
 M = 4;
 rolloff = 0.5;
 
 Rsym = 2.5e10; % symbol rate (sym/sec)
 
 span = 6; % filter span
-sps = 4; % samples per symbol
+sps = 8; % samples per symbol
 
 fs = Rsym * sps; % sampling freq (Hz)
 Tsamp = 1 / fs;
@@ -48,15 +48,20 @@ linewidthTx = 0; % Hz
 linewidthLO = 5e6; % Hz
 %linewidthLO = Rsym * 1e-3;
 
-iterations = 25;
+iterations = 1;
 avgSa = 40;
 
+TsampOrig = Tsamp;
+
 for it = 1 : iterations
   [xPSKpn, pTxLoPSK] = phaseNoise(xPSK, linewidthTx, linewidthLO, Tsamp);
   [xDEPSKpn, pTxLoDEPSK] = phaseNoise(xDEPSK, linewidthTx, linewidthLO, Tsamp);
   [xDPSKpn, pTxLoDPSK] = phaseNoise(xDPSK, linewidthTx, linewidthLO, Tsamp);
 
   for i = 1:plotlen
+    Tsamp = TsampOrig;
+    sps = 8;
+
     snr = EbN0_db(i) + 10 * log10(log2(M)) - 10 * log10(sps);
     noiseEnergy = 10 ^ (-snr / 10);
 
@@ -68,9 +73,12 @@ for it = 1 : iterations
     rDEPSK = rxFilter(yDEPSK, rolloff, span, sps);
     rDPSK = rxFilter(yDPSK, rolloff, span, sps);
 
-    rPSKSamp = rPSK(sps*span/2+1:sps:(numSymbs+span/2)*sps);
-    rDEPSKSamp = rDEPSK(sps*span/2+1:sps:(numSymbs+span/2)*sps);
-    rDPSKSamp = rDPSK(sps*span/2+1:sps:(numSymbs+span/2)*sps);
+    sps = 2;
+    Tsamp = TsampOrig * 4;
+
+    rPSKSamp = rPSK(1:2:end);
+    rDEPSKSamp = rDEPSK(1:2:end);
+    rDPSKSamp = rDPSK(1:2:end);
 
     [rPSKSampEq, phiestsPSK] = phaseNoiseCorr(rPSKSamp, M, pi/M, avgSa);
     [rDEPSKSampEq, phiestsDEPSK] = phaseNoiseCorr(rDEPSKSamp, M, 0, avgSa);
@@ -99,7 +107,7 @@ for it = 1 : iterations
 
     if EbN0_db(i) == 8 && it == 1
       figure(1234);
-      plot(repelem(-phiestsPSK, sps));
+      plot(repelem(-phiestsPSK, 8));
       hold on;
       plot(pTxLoPSK);
       legend('estimate', 'actual');
diff --git a/rxFilter.m b/rxFilter.m
index eb25a32..e8e572a 100644
--- a/rxFilter.m
+++ b/rxFilter.m
@@ -8,13 +8,23 @@ function r = rxFilter(y, rolloff, span, sps)
   %% Output:
   %%  - r: filtered signal (energy not normalized),
   %%       normalize and then sample.
-  filter = comm.RaisedCosineReceiveFilter...
+  rxfilter = comm.RaisedCosineReceiveFilter...
                  ('Shape', 'Square root', ...
                   'RolloffFactor', rolloff, ...
                   'FilterSpanInSymbols', span, ...
                   'InputSamplesPerSymbol', sps, ...
-                  'DecimationFactor', 1);
+                  'DecimationFactor', 4, ...
+                  'DecimationOffset', mod(span+1, 4));
 
-  r = filter([y; zeros(span * sps, 1)]);
-  r = r(span * sps / 2 + 1 : end); % truncate filter transients
+  coef = coeffs(rxfilter);
+  filter_fft = fft(coef.Numerator, length(y));
+  y_fft = fft(y);
+
+  rs = ifft(y_fft .* filter_fft.');
+
+  %% re-order signal due to circular conv
+  l = (length(coef.Numerator) - 1) / 2;
+  rr = [rs(l:end); rs(1:l-1)];
+
+  r = downsample(rr, 4, 3);
 end
diff --git a/splitstepfourier.m b/splitstepfourier.m
index f71828d..f0ec5fb 100644
--- 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.
diff --git a/ssf_pdm.m b/ssf_pdm.m
new file mode 100644
index 0000000..55511fe
--- /dev/null
+++ b/ssf_pdm.m
@@ -0,0 +1,61 @@
+function [xCDKerr, yCDKerr] = ...
+         ssf_pdm(x, y, 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))
+  %%  - lambda: wavelength (nm)
+  %%  - z: length of fibre (km)
+  %%  - Tsamp: sampling time (s)
+  %% Output:
+  %%  - xCDKerr: x after being dispersed.
+
+  %% 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
+
+  gamma = gamma * 1e-3; % watt^-1 / m
+  dz = 1; % km
+  dz = dz * 1e3; % m
+  stepnum = z / dz;
+
+  alpha = 0.2; % dB/km
+  alpha = alpha / 10 * log(10); % /km
+  alpha = alpha * 1e-3; % /m
+
+  hhz = 1j * (8/9) * gamma * dz;
+
+  P = abs(x) .^ 2 + abs(y) .^ 2;
+
+  xCDKerr = x .* exp(P .* hhz / 2 - alpha * dz / 4);
+  yCDKerr = y .* exp(P .* hhz / 2 - alpha * dz / 4);
+
+  for i = 1 : stepnum
+    xDFT = fft(xCDKerr);
+    yDFT = fft(yCDKerr);
+    n = length(xCDKerr);
+    fs = 1 / Tsamp;
+
+    omega = (2*pi * fs / n * [(0 : floor((n-1)/2)), (-ceil((n-1)/2) : -1)]).';
+    dispDFT = exp(-1j * omega.^2 * D * lambda^2 * dz / (4 * pi * c));
+
+    xCDKerr = ifft(xDFT .* dispDFT);
+    yCDKerr = ifft(yDFT .* dispDFT);
+
+    P = abs(xCDKerr) .^ 2 + abs(yCDKerr) .^ 2;
+    xCDKerr = xCDKerr .* exp(P .* hhz - alpha * dz / 2);
+    yCDKerr = yCDKerr .* exp(P .* hhz - alpha * dz / 2);
+    if mod(i, 50) == 0
+      xCDKerr = xCDKerr * 10; % amplification
+      yCDKerr = yCDKerr * 10;
+    end
+  end
+  P = abs(xCDKerr) .^ 2 + abs(yCDKerr) .^ 2;
+  xCDKerr = xCDKerr .* exp(-P .* hhz / 2 + alpha * dz / 4);
+  yCDKerr = yCDKerr .* exp(-P .* hhz / 2 + alpha * dz / 4);
+end
+%% References
+%% [1]: S.J. Savory, Digital filters for coherent optical receivers, 2008.
diff --git a/txFilter.m b/txFilter.m
index 507ffe1..10c387d 100644
--- a/txFilter.m
+++ b/txFilter.m
@@ -8,11 +8,17 @@ function x = txFilter(modData, rolloff, span, sps)
   %% Output:
   %%  - x: pulse-shaped waveform
 
-  filter = comm.RaisedCosineTransmitFilter...
+  txfilter = comm.RaisedCosineTransmitFilter...
                ('Shape', 'Square root', ...
                 'RolloffFactor', rolloff, ...
                 'FilterSpanInSymbols', span, ...
                 'OutputSamplesPerSymbol', sps, ...
                 'Gain', sqrt(sps)); % so that output has energy 1
-  x = filter([modData; zeros(span, 1)]);
+
+  coef = coeffs(txfilter);
+  filter_fft = fft(coef.Numerator, length(modData) * sps);
+  modData_fft = fft(upsample(modData, sps));
+  x = ifft(modData_fft .* filter_fft.');
+  l = (length(coef.Numerator) - 1) / 2;
+  x = [x(l:end); x(1:l-1)];
 end