Added adaptive CMA equalizer, "fixed" CD and phase noise

[4yp.git] / adaptiveCMA.m

diff --git a/adaptiveCMA.m b/adaptiveCMA.m
new file mode 100644 (file)
index 0000000..084eab4
--- /dev/null
+++ b/adaptiveCMA.m
@@ -0,0 +1,57 @@
+function adaptFilterOut = adaptiveCMA(rSampled)
+  %% adaptive filter
+  %% CMA
+  taps = 19; % ODD taps
+  hxx = 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;
+
+  mu = 1e-3;
+  numSymbs = length(rSampled);
+
+  %% Check average energy of symbols
+  rSampledUnitEnergy = normalizeEnergy(rSampled, numSymbs, 1);
+
+  adaptFilterOut = zeros(numSymbs, 1);
+  converged = 0;
+  convergeCount = 0;
+
+  for it = 1:numSymbs
+    if it <= (taps - 1) / 2;
+      xp = [zeros((taps - 1) / 2 - it + 1, 1); rSampledUnitEnergy(1:it + (taps - 1) / 2)];
+    elseif it + (taps - 1) / 2 > numSymbs
+      xp = [rSampledUnitEnergy(it - (taps - 1) / 2 : end); zeros(it + (taps - 1) / 2 - numSymbs, 1)];
+    else
+      xp = rSampledUnitEnergy(it - (taps - 1) / 2 : it + (taps - 1) / 2);
+    end
+
+    xout = sum(hxx .* xp);
+    ex = 1 - abs(xout) ^ 2;
+
+    if abs(ex) < 1e-3
+      convergeCount = convergeCount + 1;
+    else
+      convergeCount = 0;
+    end
+    if ~converged && convergeCount >= 10
+      converged = 1
+      it
+    end
+
+    adaptFilterOut(it) = xout;
+
+    hxx = hxx + mu * ex * xout * conj(xp);
+  end
+
+%{
+  %% try MATLAB builtin equalizer
+  alg = cma(mu);
+  eqObj = lineareq(taps, alg);
+  eqObj.Weights((taps + 1) / 2) = 1;
+  rPadded = [rSampledUnitEnergy; zeros((taps - 1) / 2, 1)];
+  matlabEq = equalize(eqObj, rPadded);
+  matlabEq = matlabEq((taps + 1) / 2 : end);
+%}
+end