Working Kerr effect; PDM; speedups; removed unused files

[4yp.git] / pdm_adaptiveCMA.m

diff --git a/pdm_adaptiveCMA.m b/pdm_adaptiveCMA.m
new file mode 100644 (file)
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