Commit | Line | Data |
---|---|---|
5117fc58 AIL |
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 | ||
79f1c8f1 | 6 | taps = 15; % Number of taps. Should be odd. |
5117fc58 AIL |
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 | |
79f1c8f1 AIL |
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)]; | |
5117fc58 AIL |
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 |