pdm_adaptiveCMA.m

   1 function [x, y] = pdm_adaptiveCMA(rx, ry)
   2   taps = 19; % ODD taps
   3   hxx = zeros(taps, 1);
   4   hxy = zeros(taps, 1);
   5   hyx = zeros(taps, 1);
   6   hyy = zeros(taps, 1);
   7   %% hxx: real indices  -K, ..., 0, ..., K.   K = floor(taps/2)
   8   %%    MATLAB indices   1      1+K     taps
   9   %% initialize hxx, hxx[0] = 1, hxx[k] = hxx[-k] = 0
  10   hxx(ceil(taps/2)) = 1;
  11   hxy(ceil(taps/2)) = 1;
  12   hyx(ceil(taps/2)) = 1;
  13   hyy(ceil(taps/2)) = 1;
  14
  15   mu = 1e-3;
  16   numSymbs = length(rx);
  17
  18   %% Check average energy of symbols
  19   rx = rx / sqrt(mean(abs(rx) .^ 2));
  20   ry = ry / sqrt(mean(abs(ry) .^ 2));
  21
  22   x = zeros(numSymbs, 1);
  23   y = zeros(numSymbs, 1);
  24   converged = 0;
  25   convergeCount = 0;
  26   convergeIdx = Inf;
  27
  28   for it = 1:numSymbs
  29     if it <= (taps - 1) / 2;
  30       xp = [zeros((taps - 1) / 2 - it + 1, 1); rx(1:it + (taps - 1) / 2)];
  31       yp = [zeros((taps - 1) / 2 - it + 1, 1); ry(1:it + (taps - 1) / 2)];
  32     elseif it + (taps - 1) / 2 > numSymbs
  33       xp = [rx(it - (taps - 1) / 2 : end); zeros(it + (taps - 1) / 2 - numSymbs, 1)];
  34       yp = [ry(it - (taps - 1) / 2 : end); zeros(it + (taps - 1) / 2 - numSymbs, 1)];
  35     else
  36       xp = rx(it - (taps - 1) / 2 : it + (taps - 1) / 2);
  37       yp = ry(it - (taps - 1) / 2 : it + (taps - 1) / 2);
  38     end
  39
  40     xout = sum(hxx .* xp) + sum(hxy .* yp);
  41     yout = sum(hyx .* xp) + sum(hyy .* yp);
  42     ex = 1 - abs(xout) ^ 2;
  43     ey = 1 - abs(yout) ^ 2;
  44
  45     if abs(ex) < 0.01
  46       convergeCount = convergeCount + 1;
  47     else
  48       convergeCount = 0;
  49     end
  50     if ~converged && convergeCount >= 10
  51       converged = 1;
  52       convergeIdx = it;
  53     end
  54
  55     x(it) = xout;
  56     y(it) = yout;
  57
  58     hxx = hxx + mu * ex * xout * conj(xp);
  59     hxy = hxy + mu * ex * xout * conj(yp);
  60     hyx = hyx + mu * ey * yout * conj(xp);
  61     hyy = hyy + mu * ey * yout * conj(yp);
  62
  63
  64     if sum(abs(hxx - hyx)) < 0.01 && sum(abs(hxy - hyy)) < 0.01
  65       hxx = 0.5 * (hxx + flipud(conj(hyy)));
  66       hyy = conj(flipud(hxx));
  67       hxy = 0.5 * (hxy - conj(flipud(hyx)));
  68       hyx = -conj(flipud(hxy));
  69     end
  70
  71   end
  72 end