Почти самый простой MIMO канал с замираниями (модель Кронекера прилагается)

Всем доброго времени суток!

Рассматривая тему пространственного разнесения ^[1], вскользь мы уже коснулись и вопросов замираний в каналах связи, и того, почему такие замирания возникают. Сегодня предлагаю поговорить об этой теме чуть более подробно.

Более того, в сегодняшней статье мы слегка коснемся вопросов пространственной корреляции MIMO каналов, которые традиционно разбираются в рамках курса профессора Хаардта ^[2] "Mobile communications" (CSP ^[3]), и поэтому, я думаю, такая публикация тем более не будет лишней. И, отмечу сразу, круг вопросов не ограничивается применением только в мобильной связи.

В общем, всех интересующихся тематикой беспроводной связи приглашаю к прочтению, и поехали!

Colin Slater "Fading memories" ^[4]. Да, сводить всё к одной теме — это, наверное, перебор, однако какой-то такой метафорой, я думаю, можно представить, что именно получает приемник после всевозможных воздействий на изначальный, информационный сигнал...

О чем вообще идет речь?

Начнем с того, что в природе существуют, так называемые, мультипликативные помехи, влияющие на принимаемую мощность сигнала — замирания (fading).

Замирания бывают быстрыми и медленными (fast and slow fading).

Рис. 1. Колебания мощности сигнала в беспроводных каналах в зависимости от расстояния. Средний уровень потерь распространения монотонно увеличивается с увеличением дальности. Локальные отклонения могут возникать из-за макроскопических (медленных) и микроскопических (быстрых) замираний [1, с.14 ].

Сознаюсь сразу, сегодня с медленными замираниями мы работать не будем, а вот про быстрые поговорим достаточно подробно.

Рис. 2. Иллюстрация эффекта затенения: ^[5] как правило, медленные замирания возникают из-за затенений (shadowing). Для моделирования обыкновенно выбирают логнормальное распределение ^[6].

Желающих почитать немного про медленные замирания, в частности, и про то, как можно моделировать спутниковые каналы, в целом, приглашаю сюда ^[7]!

Быстрые замирания

Быстрые замирания возникают, как правило, по двум основным причинам:

• из-за уже упомянутого нами многолучевого распространения (multipath propagation) и/или
• из-за Допплеровских сдвигов частоты.

Последние очень заметны для спутниковых систем ^[8], однако не всегда учитываются для случаев наземной связи из-за того, что базовые станции в данном случае, как правило, неподвижны.

Рис. 3. Иллюстрация причин возникновения замираний. ^[9] Отдельное внимание уделено рассеиванию (scattering) и дифракции: именно в силу данных явлений и формируется, так называемое, многолучевое распространение, приводящее к межсимвольным помехам (ISI ^[10]).

Но и это ещё далеко не всё.

Selective fading vs. Flat fading

Выше мы разделили наши помехи по характеру возникновения. Однако, помехи можно разделить ещё и по характеру воздействия на передаваемый сигнал. И здесь нам понадобится понятие избирательности канала.

Приведем небольшую классификацию по [1, c. 14-18]. Итак, быстрые замирания могут быть:

1. Избирательными (selective)
   а. Частотно избирательными (frequency selective)
   б. Избирательными во временной области (time selective)
   в. Пространственно избирательными (это относится к вопросу об углах прихода и отправки ЭМ волн — сегодня мы этот вопрос разбирать не будем)
2. Плоскими (flat) — тяготеющими больше к характеру медленных замираний (да, вот такой вот парадокс)

Что подразумевает последний термин, объясним от обратного.

Во-первых, если ширина спектра передаваемого нами сигнала превышает когерентную полосу пропускания канала ^[11], то говорят, что мы имеем дело с частотно-селективными помехами.

Рис. 4. Иллюстрация когерентной полосы пропускания. Думаю, не будет моветоном прикрепить иллюстраций из Вики ^[12], если сам их туда и добавил.

Обратите внимание на переменную Delay spread — разброс задержек. Именно этот разброс в задержках между приходом разных копий одного сигнала и измеряют, когда определяют характеристики того или иного реального канала.

Рис. 5. Типичный профиль задержки (мощности) — средняя мощность как функция задержки .[1, c. 16]

Рис. 6. Частотно-селективный, но не избирательно по времени канал ^[13].

Во-вторых, если длина нашего сигнала во временной области превышает время когерентности канала ^[17], то говорят, что мы имеем дело с каналом, изменяющимся во времени.

Рис. 7. Иллюстрация времени когерентности. Обратите внимание, здесь максимальная допплеровская частота отражает движение самого мобильного терминала.

Рис. 8. Почти неизбирательный по частоте, но избирательный по времени: быстро затухающий канал. ^[13]

Ну, и соответственно, если нам удастся каким-то чудом избежать вышеперечисленного, то мы придем к самомому простому и удобному случаю — к плоским замираниям [1, c. 14-18][2, c. 88-92].

Самый простой канал с замираниями

Возможно, вы где-то уже встречали подобный способ моделирования плоского канала с замираниями:

h = (1/sqrt(2))*(randn(size(msg)) + 1j*randn(size(msg))) % MatLab / Octave

Что подразумевается:

• отсутствие прямой видимости (вполне оправданный для наземной связи случай)
• узкополосный канал (narrowband channel) — вписались в когерентную полосу
  • кстати, пример широкополосного канала мы рассматривали в данной статье (рис. 4). ^[19]
• частотной селективности тоже нет — вписались в когерентное время
• отсутствие затенений и Допплеровских сдвигов частоты

То есть вот он собственной персоной: плоский канал с замираниями. Для случая без прямой видимости канал ещё называют к тому же рэлеевским.

Лорд Рэлей

Из чего проистекает такая модель?

Во-первых, говоря о современных системах связи, мы скорее всего подразумеваем работу с цифровыми сигналами. Более того, работаем мы к тому же скорее всего c квадратурными типами модуляции (I/Q ^[20]), то есть с сигналами комплексными.

Рис. 9. Сигнальное созвездие QPSK ^[21]. Полярные координаты, единичная окружность. Например, символу 11 будет соответствовать координата равная примерно 0,7+0,7i.

Логично предположить воздействие помех как на реальную, так и на мнимую составляющую. Более того, стоит вслед за Центральной предельной теоремой ^[22] предположить также гауссовский характер таких воздействий. Запишем формулой [1, с. 39]:

где и — это нормально распределенные случайные величины.

Чаще всего используется модель с нормированной средней мощностью:

А значит .

Если мы предполагаем независимые искажения составляющих огибающей сигнала: как синфазной (I — In-phase), так и квадратурной (Q — Quadrature), — то положение числа, на которое будет перемножен изначальный символ (отсюда мультипликативность помехи), в полярных координатах может быть вычислено буквально по теореме Пифагора:

где и — это тоже нормально распределенные случайные величины. Для нормированной мощности , .

И вот перед нами формула рэлеевского процесса ^[23]! [2, с.78 ]

Рис. 10. Гауссовские генераторы в квадратуре для моделирования замираний Рэлея и Райса [3, p.125]. О райсовском канале мы поговорим чуть позже.

Если мы рассматриваем самый простой случай без какой-либо пространственной корреляции между каналами (допустим, что у нас есть несколько пространственных каналов — MIMO), то говорят, что канал пространственно белый — spatially white. А так как все каналы независимы, то, следовательно, и "генерить" можно не только случайные величины, распределенные по Рэлею, но и вектора, и матрицы, и тензоры [3, c.125 ]:

где и — это матрицы, состоящие из нормально распределенных значений.

Откуда взялась 1/2 под корнем?

Тут всё довольно просто: как мы уже сказали выше, мощность отклика канала должна быть равна единице, т.е.:

Именно за этим и вводят масштабирующий фактор

Продемонстрируем (возьмем python для разнообразия):

import numpy as np

N = int(1e6)
h = 1/np.sqrt(2)*(np.random.randn(N) + 1j*np.random.randn(N))
P = np.mean(h*np.conj(h))
print(np.round(P))

>>> (1+0j)

Подытожим:
Помеховая компонента при отсутствии прямой видимости может быть описана как n-мерный вектор (в зависимости от размерности канала) независимых идентично распределенных (IID — independent identically distributed) по комплексному гауссовскому закону с круговой симметрией и нулевым мат. ожиданием (ZMCSCG — zero-mean circularly symmetric complex Gaussian) чисел.

Целая скороговорка!

Почти самый простой канал с замираниями

Хорошо, а если прямая видимость всё-таки есть? Может быть, есть на этот счет уже готовые математические модели?

Конечно, есть! Например такая [4] ^[25], обобщенная для случая с несколькими приёмными и передающими антеннами (MIMO):

где — это собственно и есть комплексная огибающая импульсной характеристики канала (или канальная матрица, если речь идет о MIMO ^[26]), — это компонента прямой видимости (Line-of-Sight), а — это компонента отсутствия прямой видимости (Non-Line-of-Sight), её мы уже рассмотрели.

Обратите внимание на переменную K — это так называемый райсовский фактор (Rician/Ricean factor) — отношение мощности прямого пути распространения к мощности всех остальных путей. Именно он определяет, насколько сильно будет подвержен наш канал помехам.

Stephen O. "Steve" Rice ^[27]

Слагаемое представляет из себя математическое ожидание канальной матрицы и моделироваться может исходя из геометрии приемных и передающих массивов антенн:

где и — это отклики приемного и передающего массивов антенн сответственно, а и — это углы прихода (angle of arrival) и отправки (angle of departure).

Множители, в общем-то, симметричные и могут быть описаны общей формулой:

где — это расстояние между антеннами (antenna spacing) в длинах волн, а — это количество элементов в массиве антенн (рис. 11).

Рис. 11. Геометрия линейного массива ^[28]. Для случая Massive MIMO геометрия, конечно, будет уже другой. Подробнее об этом прочитать можно, например, здесь ^[29].

Ну, а теперь, когда модель хоть сколько-то разобрана, можно её и упростить.

Для случая SISO (самого простого и базового) формула (5) примет вид:

ибо для 1 передающей и 1 приемной антенны .

И это, как нельзя кстати, подходит под модель на рисунке 12.

Рис. 12. Схема формирования канала Райса (узкополосный сигнал) [3, p.127]. При канал будет всё больше вырождаться в Рэлеевский. Поэтому я и употребил слово почти в подзаголовке: канал Рэлея — это частный случай канала Райса. Это также видно из формулы (5).

Поясним: — это райсовский масштабирующий параметр (Rician scale parameter), а — это параметр не-центральности (Rician noncentrality parameter).

Моделирование почти самого простого канала с замираниями в MatLab

Сверим полученные выкладки с чем-нибудь уже проверенным. Например, со специальными функциями MatLab ^[32].

Для проверки модели со своей стороны сгенерируем случайное двоичное сообщение (длина сообщения равна 100000 бит ), промодулируем его с помощью M-PSK ^[21] или M-QAM ^[33], поэлементно умножим с процессом замираний, добавим белый гауссовский шум, выровняем эквалайзером Zero-Forcing, демодулируем и рассчитаем битовую ошибку (BER — bit error ratio) (рис. 13). Количество испытаний будет равным ста.

Рис.13. Блок-схема модели системы передачи, которую опишет наш скрипт.

clear all; close all; clc 

EbNo = 0:40; 
K = [4.0; 0.6]; 
M = [4; 8; 16; 64; 256]; %Positions of modulation (M-PSK or M-QAM) 

for k = 1:length(K) 
    for m = 1:length(M) 
        message = randi([0, M(m)-1], 100000, 1); 
        if M(m) >= 16 
            mod_msg = qammod(message, M(m), pi/4, 'gray'); 
            ric_ber(:, m, k) = berfading(EbNo,'qam',M(m),1,K(k)); 
        else 
            mod_msg = pskmod(message, M(m), pi/4, 'gray'); 
            ric_ber(:, m, k) = berfading(EbNo, 'psk', M(m), 1, K(k)); 
        end 
        Es = mean(abs(mod_msg).^2); 
        No = Es./((10.^(EbNo./10))*log2(M(m))); 

        h = sqrt( K(k)/(K(k)+1)) +... 
        sqrt( 1/(K(k)+1))*(1/sqrt(2))*(randn(size(mod_msg))...
         + 1j*randn(size(mod_msg))); 
        ric_msg = mod_msg.*h; % Rician flat fading 

        for c = 1:100 
            for jj = 1:length(EbNo) 
                noisy_mod = ric_msg +... 
                sqrt(No(jj)/2)*(randn(size(mod_msg))+...
                1j*randn(size(mod_msg))); %AWGN 
                noisy_mod = noisy_mod ./ h; % zero-forcing equalization 
                if M(m) >= 16 
                    demod_msg = qamdemod(noisy_mod, M(m), pi/4, 'gray'); 
                else 
                    demod_msg = pskdemod(noisy_mod, M(m), pi/4, 'gray'); 
                end 
                [number,BER(c,jj)] = biterr(message,demod_msg); 
            end 
        end 
        sum_BER(:,m, k) = sum(BER)./c; 
    end 
end 

figure(1) 

semilogy(EbNo, sum_BER(:,1,1), 'b-o', EbNo, sum_BER(:,2,1), 'r-o',... 
EbNo, sum_BER(:,3,1), 'g-o', EbNo, sum_BER(:,4,1), 'c-o',...
EbNo, sum_BER(:,5,1), 'k-o',... 
EbNo, ric_ber(:,1,1), 'b-', EbNo, ric_ber(:,2,1), 'r-',...
EbNo, ric_ber(:,3,1), 'g-', EbNo, ric_ber(:,4,1), 'c-',...
EbNo, ric_ber(:,5,1), 'k-', 'LineWidth', 1.5) 
title('Rician model (K = 4.0)') 
legend('QPSK(simulated)', '8-PSK(simulated)',... 
'16-QAM(simulated)', '64-QAM(simulated)' ,'256-QAM(simulated)',... 
'QPSK(theory)','8-PSK(theory)', '16-QAM(theory)',... 
'64-QAM(theory)' ,'256-QAM(theory)','location','best') 
xlabel('EbNo (dB)') 
ylabel('BER') 
grid on 

figure(2) 

semilogy(EbNo, sum_BER(:,1,2), 'b-o', EbNo, sum_BER(:,2,2), 'r-o',... 
EbNo, sum_BER(:,3,2), 'g-o', EbNo, sum_BER(:,4,2), 'c-o',...
EbNo, sum_BER(:,5,2), 'k-o',... 
EbNo, ric_ber(:,1,2), 'b-', EbNo, ric_ber(:,2,2), 'r-',...
EbNo, ric_ber(:,3,2), 'g-', EbNo,ric_ber(:,4,2), 'c-',...
EbNo, ric_ber(:,5,2), 'k-','LineWidth', 1.5)
title('Rician model (K = 0.6)') 
legend('QPSK(simulated)', '8-PSK(simulated)',... 
'16-QAM(simulated)', '64-QAM(simulated)' ,'256-QAM(simulated)',... 
'QPSK(theory)','8-PSK(theory)',... 
'16-QAM(theory)', '64-QAM(theory)' ,'256-QAM(theory)','location','best') 
xlabel('EbNo (dB)') 
ylabel('BER') 
grid on

Должно получится что-то такое:

Рис. 14. Коэффициент битовых ошибок для рассмотренной модели (К = 4.0).

Fig. 15. Коэффициент битовых ошибок для рассмотренной модели (K = 0.6).

Не без огрехов, конечно (из-за ограниченности по количеству усреднений), но в целом, видно, что работает!

Почти самый простой канал с замираниями и пространственная корреляция

А теперь, завершая наше путешествие в мир статистических выкладок, давайте немного приблизимся к реальности.

Выше мы рассматривали случай, когда компонента отсутствия прямой видимости была пространственно белой. Однако, от реальных систем стоит всё-таки ждать некоторой корреляции между антеннами, а значит и между путями распространения ЭМ волны.

Введем новую переменную — матрицу пространственной корреляции:

где vec обозначает операцию векторизации ^[36]. Влиять эта матрица будет, кстати, только на компоненту отсутствия прямой видимости ^[37], и значит в данном случае .

Резонно предположить, что в реальной жизни сначала через замеры определяют данную матрицу, а уже потом используют её для моделирования:

где — это IID матрица (w означает spatially white).

Обратите внимание, возведение матрицы в степень — это не возведение в степень каждого её элемента. Следовательно в нашем случае в MatLab нельзя использовать функцию sqrt(). Использовать нужно оператор ^.
В случае Python нельзя использовать оператор возведения в степень ** и метод sqrt популярного модуля numpy. Использовать нужно метод fractional_matrix_power ^[38] модуля scipy.

Если приёмник и передатчик достаточно независимы друг от друга, как на рисунке 16, то можно применить модель Кронекера [1, c. 40], чтобы уменьшить количество коэффициентов для расчета:

где обозначает произведение Кронекера ^[39].

Леопольд Кронекер

Отдельно рассмотрим составляющие выражение матрицы:

• — это матрица корреляции приемных антенн,

• — это матрица корреляции передающих антенн

и — это количество приемных и передающих антенн соответственно. Обе матрицы являются положительными полуопределенными ^[40] (positive semi-definite) Эрмитовыми ^[41] (Hermitian) матрицами. Более того, ранг матрицы ограничен значением .

Рис. 16.
Модель рассеивателей в случае некоррелирующих приёмника и передатчика. Зеленые точки означают рассеиватели, синие и красные — приемник и передатчик. Прошу простить за качество, это скриншот из одной из домашек из Ильменау

Канальную матрицу в данном случае можно записать как:

Однако, для ситуаций, как, например, на рисунке 17, использовать придется всё же полную корреляционную матрицу (см. формулу (9)).

Рис. 17. Модель рассеивателей в случае коррелирующих приёмника и передатчика. Зеленые точки означают рассеиватели, синие и красные — приемник и передатчик.

Дабы закрепить узнанное и посмотреть, насколько осложнит жизнь пространственная корреляция, давайте попробуем посчитать пропускную способность каналов. Например, для MIMO 2x2.

Для этого возьмем следующие случаи (имеем в виду всё те же плоские замирания, для простоты также рассмотрим только рэлеевский канал):

1) Пространственно белый канал;
2) Модель Кронекера с относительно слабой корреляцией:

3) Модель Кронекера с относительно сильной корреляцией:

Использовать в данной части будем Python 3.4. Функцию для расчета пропускной способности возьмем из одной из моих прошлых публикаций ^[26].

import numpy as np
from numpy import linalg as LA
import matplotlib.pyplot as plt
from scipy.linalg import fractional_matrix_power

def openloop_capacity(H_chan, SNR_dB):
    SNR = 10**(SNR_dB/10)
    Mt = np.shape(H_chan)[1]
    H_sq = np.dot(H_chan,np.matrix(H_chan, dtype=complex).H)
    lambdas = LA.eigvals(H_sq) 
    lambdas = np.sort(lambdas)[::-1]
    c = 0
    for eig in lambdas:
        c = c + np.log2(1 + SNR*eig/Mt)
    return np.real(c)

Определим основные параметры системы:

Mr = 2 # number of receive antennas
Mt = 2 # number of transmit antennas

Rr1 = np.array([[1., .4], [.4, 1.]]) # case 1
Rr2 = np.array([[1., .6], [.6, 1.]]) # case 2
Rt1 = np.array([[1., .4], [.4, 1.]]) # case 1
Rt2 = np.array([[1., .6], [.6, 1.]]) # case 2

Rr1_half = fractional_matrix_power(Rr1, 0.5)
Rr2_half = fractional_matrix_power(Rr2, 0.5)

Rt1_half = fractional_matrix_power(Rt1, 0.5)
Rt1_half_H = np.conj(Rt1_half).T

Rt2_half = fractional_matrix_power(Rt2, 0.5)
Rt2_half_H = np.conj(Rt2_half).T

Запускаем моделирование:

counter = 1000
SNR_dBs = [i for i in range(1, 21)]

C_MIMO_uncorr = np.empty((len(SNR_dBs), counter))
C_MIMO_kron1 = np.empty((len(SNR_dBs), counter))
C_MIMO_kron2 = np.empty((len(SNR_dBs), counter))
C_MIMO_kron3 = np.empty((len(SNR_dBs), counter))
C_MIMO_kron4 = np.empty((len(SNR_dBs), counter))

for c in range(counter):
    H_uncorr = (np.random.randn(Mr,Mt) + 1j*np.random.randn(Mr, Mt))/np.sqrt(2)
    H_kron1 = np.dot(Rr1_half, np.dot(H_uncorr, Rt1_half_H))
    H_kron2 = np.dot(Rr2_half, np.dot(H_uncorr, Rt2_half_H))

    for idx, SNR_dB in enumerate(SNR_dBs):
        C_MIMO_uncorr[idx, c] = openloop_capacity(H_uncorr, SNR_dB)
        C_MIMO_kron1[idx, c] = openloop_capacity(H_kron1, SNR_dB)
        C_MIMO_kron2[idx, c] = openloop_capacity(H_kron2, SNR_dB)

C_MIMO_uncorr_erg = np.mean(C_MIMO_uncorr, axis=1)
C_MIMO_kron1_erg = np.mean(C_MIMO_kron1, axis=1)
C_MIMO_kron2_erg = np.mean(C_MIMO_kron2, axis=1)

plt.figure(figsize=(7, 5), dpi=100)

plt.plot(SNR_dBs, C_MIMO_uncorr_erg,'g-o', label='Uncorrelated')
plt.plot(SNR_dBs, C_MIMO_kron1_erg, 'm-o', label='Kronecker model, case 1')
plt.plot(SNR_dBs, C_MIMO_kron2_erg, 'r-o', label='Kronecker model, case 2')

plt.title("Rayleigh flat fading (MIMO 2x2)")
plt.xlabel('SNR (dB)')
plt.ylabel('Ergodic capacity (bps/Hz)')
plt.legend()
plt.minorticks_on()
plt.grid(which='major')
plt.grid(which='minor', linestyle=':')
plt.show()

Получаем:

Рис. 18. Кривые пропускной способности для указанных экспериментов.

Что ж: разница хоть и не критическая, но есть!

Послесловие

Возможно, кто-то из читателей мог заметить, что почти всё, что мы рассмотрели сегодня (подробно и вскользь) весьма неплохо ложится на параметры матлабовского объекта comm.MIMOChannel ^[42]. Мне этот факт кажется очень важным:

• во-первых, это дополнительно систематизирует выше описанное,
• а, во-вторых, подчеркивает то, что вопросы мы рассматриваем актуальные.

Ну, и теперь, если вы будете использовать подобные блоки, то, надеюсь, будет более понятным зачем и какие ручки нужно крутить.

Всех благодарю за внимание!

P.S.

Добавлю сюда ссылку на слайды ^[43] по теме от Reiner S. Thomä ^[44] (Ильменау), по которым он нам читал в своё время лекции. Думаю, кому-нибудь пригодится!

Литература

1. Paulraj, Arogyaswami, Rohit Nabar, and Dhananjay Gore. Introduction to space-time wireless communications. Cambridge university press, 2003.
2. Goldsmith A. Wireless communications. – Cambridge university press, 2005.
3. Fontæn, Fernando Pærez, and Perfecto Mariæo Espiæeira. Modelling the wireless propagation channel: a simulation approach with Matlab. Vol. 5. John Wiley & Sons, 2008.
4. Farrokhi, Farrokh R., et al. "Spectral efficiency of FDMA/TDMA wireless systems with transmit and receive antenna arrays." IEEE transactions on wireless communications 1.4 (2002): 591-599.

Автор: Владимир Фадеев

Источник ^[45]