pdm_adaptiveCMA.m

   1 function [x, y] = pdm_adaptiveCMA(rx, ry)
   2   %% Perform adaptive equalization using CMA.
   3   %% Input: rx, ry: Both polarizations of received signal
   4   %% Output: x, y: Equalizaed signal
   5
   6   taps = 15; % Number of taps. Should be odd.
   7   mu = 1e-3; % Convergence parameter for gradient descent.
   8
   9   hxx = zeros(taps, 1);
  10   hxy = zeros(taps, 1);
  11   hyx = zeros(taps, 1);
  12   hyy = zeros(taps, 1);
  13   %% hxx: real indices  -K, ..., 0, ..., K.   K = floor(taps/2)
  14   %%    MATLAB indices   1      1+K     taps
  15
  16   %% Initialize hxx, hxx[0] = 1, hxx[k] = hxx[-k] = 0
  17   hxx(ceil(taps/2)) = 1;
  18   hxy(ceil(taps/2)) = 1;
  19   hyx(ceil(taps/2)) = 1;
  20   hyy(ceil(taps/2)) = 1;
  21
  22   numSymbs = length(rx);
  23
  24   %% Normalize to unit power.
  25   rx = rx / sqrt(mean(abs(rx) .^ 2));
  26   ry = ry / sqrt(mean(abs(ry) .^ 2));
  27
  28   x = zeros(numSymbs, 1);
  29   y = zeros(numSymbs, 1);
  30
  31   %% Run CMA twice so that the first symbols were also equalized
  32   for loops = 1:2
  33     %% Loop through each symbol
  34     for it = 1:numSymbs
  35       %% Construct block of length equal to filter length (taps)
  36       if it <= (taps - 1) / 2;
  37         %% If near the start, prepend zeros
  38         xp = [zeros((taps - 1) / 2 - it + 1, 1); rx(1:it + (taps - 1) / 2)];
  39         yp = [zeros((taps - 1) / 2 - it + 1, 1); ry(1:it + (taps - 1) / 2)];
  40       elseif it + (taps - 1) / 2 > numSymbs
  41         %% If near the end, append zeros
  42         xp = [rx(it - (taps - 1) / 2 : end); ...
  43               zeros(it + (taps - 1) / 2 - numSymbs, 1)];
  44         yp = [ry(it - (taps - 1) / 2 : end); ...
  45               zeros(it + (taps - 1) / 2 - numSymbs, 1)];
  46       else
  47         %% Just slice the signal
  48         xp = rx(it - (taps - 1) / 2 : it + (taps - 1) / 2);
  49         yp = ry(it - (taps - 1) / 2 : it + (taps - 1) / 2);
  50       end
  51
  52       %% Filtering
  53       xout = sum(hxx .* xp) + sum(hxy .* yp);
  54       yout = sum(hyx .* xp) + sum(hyy .* yp);
  55       x(it) = xout;
  56       y(it) = yout;
  57
  58       %% Caculate error
  59       ex = 1 - abs(xout) ^ 2;
  60       ey = 1 - abs(yout) ^ 2;
  61
  62       %% Update filter by gradient descent
  63       hxx = hxx + mu * ex * xout * conj(xp);
  64       hxy = hxy + mu * ex * xout * conj(yp);
  65       hyx = hyx + mu * ey * yout * conj(xp);
  66       hyy = hyy + mu * ey * yout * conj(yp);
  67
  68       %% If both filters converge to the same polarization,
  69       % re-initialize the filters.
  70       if sum(abs(hxx - hyx)) < 0.01 && sum(abs(hxy - hyy)) < 0.01
  71         hxx = 0.5 * (hxx + flipud(conj(hyy)));
  72         hyy = conj(flipud(hxx));
  73         hxy = 0.5 * (hxy - conj(flipud(hyx)));
  74         hyx = -conj(flipud(hxy));
  75       end
  76     end
  77   end
  78 end