Monday, March 10, 2014

Rotated constellation demapper design

    Hi! I have been working on an OFDM demodulator. It uses rotated constellations (Q-delay, just like in DVB-T2) to improve BER performance over frequency selective channels. I am going to share my thoughts about demapping algorithm here and describe some details. You can use MATLAB code to repeat the simulation and check results by yourself.

    Here is a structure of my work:
  1. Demapping. LLR calculation.
  2. Rotated constellations. Q-delay.
  3. Rotated constellation demapper.
    1. Direct implemention
    2. QRD based demapper



        1. Demapper
    Demappers are used in receivers to decide what bit was transmitted according to received symbol (for example QAM or PSK) and know ideal signal constellation. There are two types of demappers:
  • Hard decision demapper. Decides what exact bits were transmitted.
  • Soft decision demapper. Evaluates Logarithmic Likelihood Ratio (LLR) for each bit.
    I am going to use the second one. Exact LLR for a bit b is defined as
where is a set of constellation points, that defines b = 0, and  defines b = 1; x is a received symbol; is a noise variance. To estimate approximate LLR I can use only the nearest constellation points from sets  and  [1]:
.
    Now I can remove logarithmic operations and simplify LLR calculation:
    LLR values are used as soft decoder inputs. I am going to use LDPC decoder in the latter examples.

    Here is an example of QAM 16 mapper and soft demapper implementations in MATLAB.

%% mapper
m = 16;
n = 1e4;

bitspersymbol = log2(m);

constel = qammod(0:m-1, m, 0, 'gray');
constel = constel / sqrt(mean(abs(constel).^2));

dat = randi([0 m-1], 1, n);
mappeddat = constel(dat + 1);

%% awgn channel
SNR = 17;
variance = 1/10^(SNR/20);
chandat = mappeddat + variance*(randn(1,n) + 1i*randn(1,n))/sqrt(2);

%% demapper
S0 = zeros(bitspersymbol, m/2);
S1 = zeros(bitspersymbol, m/2);

for b = 1:bitspersymbol
    point_id = bitget(0:m-1, b);
    S1(b, :) = constel(point_id == 1);
    S0(b, :) = constel(point_id == 0);
end

LLR = zeros(bitspersymbol, n);
for b = 1:bitspersymbol
    [~, s0_id] = min( abs( repmat(S0(b, :).', 1, n) - repmat(chandat, m/2, 1) ) );
    [~, s1_id] = min( abs( repmat(S1(b, :).', 1, n) - repmat(chandat, m/2, 1) ) );

    s0 = S0(b, s0_id);
    s1 = S1(b, s1_id);
   
    LLR(b, :) = 1/variance^2 * (abs(chandat - s1).^2 - abs(chandat - s0).^2);
end

hard_decision = 2.^(0:bitspersymbol-1) * (LLR<0);

%% plots
figure(1)
plot(chandat, '.')
hold on
plot(constel, '.g')
hold off
legend('Received symbols', 'Ideal constellation positions')

figure(2)
for b = 1:bitspersymbol
    subplot(bitspersymbol/2, 2, b)
    plot(LLR(b,:),'.')
    title([num2str(b) ' bit LLR'])
end





Output plots for 20 dB and 30 dB SNR.

    You can see that LLR values highly depend on SNR, and on bit position. Some of the bits are stronger than others. This fact should be taken into account when you design a system, because if you use irregular code, some of the encoded bits could be more important.

    If there is a frequency selective fading the SNR would be different for each subcarrier. Assume that we have white Gaussian noise at the receiver side. Before the equalization noise variance is constant and symbol (subcarrier) power depends on the channel condition. Equalizer estimates channel frequency response and compensates gain of the channel (equalizes subcarriers powers). Now noise variance becomes variable.


    In this case LLR on each subcarrier should be calculated using proper noise variance value. I prefer to use a demapper with variable reference points. In this case noise variance remains same for each subcarrier.


        2. Rotated constellation with Q-delay

    It is very simple and effective technique. Mapped symbol is rotated by fixed angle, which depends only on a modulation order. In the following table you can see rotation angles for different modulations of DVB-T2 standard [2].

Modulation
Rotation Angle [degree]
QPSK
29
16-QAM
16.8
64-QAM
8.6
256-QAM
atan(1/16)

%% mapper
m_id = 1;
n = 1e4;
m_table = [4 16 64 256];
rot_ang_table = [[29 16.8 8.6]/180*pi atan(1/16)];

m = m_table(m_id);
rot_ang = rot_ang_table(m_id);

bitspersymbol = log2(m);

PAM_map = pammod(0:sqrt(m)-1, sqrt(m), 0, 'gray');
QAM_map = (repmat(PAM_map.', 1, sqrt(m)) + 1i*(repmat(PAM_map, sqrt(m), 1))).';
QAM_map =  QAM_map(:).' * exp(1i*rot_ang);

dat = randi([0 m-1], 1, n);

mappeddat = QAM_map(dat + 1);


    After rotation we should perform Q-delay operation. Each constellation point is a sum of In-phase (I) and Quadrature (Q) components: S(i) = I(i) + j*Q(i). After rotation Q component is shifted so that resulting symbol becomes a sum of two independent components S(i) = I(i) + j*Q(i-1). So, the transmitted constellations become very unusual. Number of the new possible symbol positions is M2.


New constellations construction principle




Resulting constellations


mappeddat = real(mappeddat) + 1i*imag(mappeddat([end 1:end-1]));




    A bit is transmitted by two different subcarriers and if one of them suffers deep fading, another one would save the information. This is the main idea of the technique. Also, its very important feature is that it doesn’t impair BER performance of a transmission over AWGN channel.



        3. Rotated constellation demapper
        3.1. Direct implementation

    At the receiver we can perform reverse operations, which are Q-delay and rotation. If there is no frequency selective fading it wouldn’t have any performance degradation or improvement. In other case we should take into account channel gain.

%% channel
norm_pow_gain = sqrt(mean(abs(QAM_map).^2));

channel = 1/norm_pow_gain * (randn(1,n) + 1i*randn(1,n))/sqrt(2);
chandat = mappeddat.*channel;

SNR = 30;
sigma = 1/10^(SNR/20);
chandat = chandat + sigma*(randn(1,n) + 1i*randn(1,n))/sqrt(2);


    In this model I use uncorrelated Rayleigh channel. I also assume that I have perfect channel estimation on the receiver side. First of all, I compensate phase rotation caused by the channel and shift back Q component.

%% demapper
demapper_idat = chandat.*exp(-1i*angle(channel));
demapper_idat = real(demapper_idat) + 1i*imag(demapper_idat([2:end 1]));


    Then I use the channel gain to find reference constellation points for each subcarrier to calculate LLR.

LLR = zeros(bitspersymbol, n);

for b = 1:bitspersymbol
    point_id = bitget(0:m-1, b);
    S0 = QAM_map(point_id == 0);
    S1 = QAM_map(point_id == 1);
   
    S0 = repmat(S0.', 1, n);
    S1 = repmat(S1.', 1, n);
   
    S0 = real(S0).*repmat(abs(channel),m/2,1) + ...
        1i*imag(S0).*repmat(abs(channel([2:end 1])),m/2,1);
    S1 = real(S1).*repmat(abs(channel),m/2,1) + ...
        1i*imag(S1).*repmat(abs(channel([2:end 1])),m/2,1);
   
    s0_dist = min( abs(S0 - repmat(demapper_idat, m/2, 1)).^2 );
    s1_dist = min( abs(S1 - repmat(demapper_idat, m/2, 1)).^2 );

    LLR(b, :) = 1/sigma^2 * (s1_dist - s0_dist);
end


Examples of rotated 16-QAM constellations at Rx with various channel gains

        3.2. QRD based demapper

    I desribe signal using matrix notation below. QAM symbol can be presented as an array of two real numbers x = {I; Q}. The complex rotation can be considered as a multiplication of rotation matrix  and a symbol vector x.
where  is a rotation angle. Channel attenuates both of the components independently, so the channel matrix here (assuming that we have perfect channel estimation) can be considered as a diagonal matrix.
    Additive noise can be omitted for now. At the demapper we can assume that there is some H’ channel matrix which includes constellation rotation and the channel gain.
Then the received symbols would be as follows
    So, the aim here is just to solve the set of two equations. There are some methods which can be used to find x. I use here an ML method with reduced complexity based on QR-decomposition.

    At the receiver side we know the H' matrix and the y array. Now we decompose H’ matrix into a product of orthogonal matrix Q and upper triangular matrix R. 
    Then we multiply each side by a transposed matrix Q.
    Now I consider both of components as PAM symbols and I use kind of a sphere decoder with only two layers.


    For example, to find a nearest point on 64-QAM constellation I could use 8 parallel branches, which consider all possible Qx values and searches for the nearest Ix value. Search is implemented over the real (non-complex) numbers so it can be done with 3 stage comparison operations. After all 8 nearest  values are found I calculate 8 Euclidean distances (in complex signal space domain) to find the smallest one.

    To calculate LLR of each bit I divide them into 2 groups. Each of them is independently mapped onto PAM constellation. At the receiver side they are demapped separately so that I have to calculate only  Euclidean distances (M is a QAM order).

%% demapper
demapper_idat = chandat.*exp(-1i*angle(channel));
demapper_idat = real(demapper_idat) + 1i*imag(demapper_idat([2:end 1]));

LLR = zeros(bitspersymbol, n);

H = abs(channel) + 1i*abs(channel([2:end 1]));
rot = [cos(rot_ang) -sin(rot_ang); sin(rot_ang) cos(rot_ang)];

LLR_ref = zeros(n, log2(m));
for i = 1:n
    % 1st half
    H_ = [real(H(i)) 0; 0 imag(H(i))] * rot;
    [Q, R] = qr(H_);
    Y_ = Q'*[real(demapper_idat(i)); imag(demapper_idat(i))];
   
    dist1 = abs(Y_(2) - PAM_map*R(2,2)).^2;

    dist2 = zeros(1, sqrt(m));
    for j = 1:sqrt(m)
        dist2(j) = min( abs( Y_(1) - PAM_map*R(1,1) - PAM_map(j)*R(1,2) ) )^2;
    end
    dist = dist1 + dist2;
   
    for bit_id = 1:log2(sqrt(m))
        id0 = bitget((0:sqrt(m)-1), bit_id) == 0;
        id1 = bitget((0:sqrt(m)-1), bit_id) == 1;
        LLR(bit_id, i) = 1/sigma^2 * min(dist(id1)) -  min(dist(id0));
    end
   
    % 2nd half
    H_ = [imag(H(i)) 0; 0 real(H(i))] * rot';
    [Q, R] = qr(H_);
    Y_ = Q'*[imag(demapper_idat(i)); real(demapper_idat(i))];
   
    dist1 = abs(Y_(2) - PAM_map*R(2,2)).^2;

    dist2 = zeros(1, sqrt(m));
    for j = 1:sqrt(m)
        dist2(j) = min( abs( Y_(1) - PAM_map*R(1,1) - PAM_map(j)*R(1,2) ) )^2;
    end
    dist = dist1 + dist2;
   
    for bit_id = 1:log2(sqrt(m))
        id0 = bitget((0:sqrt(m)-1), bit_id) == 0;
        id1 = bitget((0:sqrt(m)-1), bit_id) == 1;
        LLR(bit_id+log2(sqrt(m)), i) = 1/sigma^2 * min(dist(id1)) - min(dist(id0));
    end
end
    The QR decomposition of the real value matrix can be implemented with a few CORDIC cores.

    The demapping algorithm is very flexible and scalable.
  • It is easy to pipeline demapper’s structure for maximal throughput;
  • It can be efficiently implemented for variable order of QAM constellations;
  • I suppose that it could be further improved for specified terms.
        References

        [1] Viterbi, A. J., "An Intuitive Justification and a Simplified Implementation of the MAP Decoder for Convolutional Codes," IEEE Journal on Selected Areas in Communications, vol. 16, No. 2, pp 260–264, Feb. 1998

        [2] DVB-T2 standard (pdf)

Creative Commons License
This work by Georgiy Trofimov is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License.
Permissions beyond the scope of this license may be available at https://www.blogger.com/profile/13302434880547879271.

4 comments: