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

[4yp.git] / adaptiveCMA.m

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