Хабрахабр

[Из песочницы] Моделируем алгоритм MUSIC для задач определения направления прихода электромагнитной волны

aaspcats

Предисловие

Давным-давно, в далеких 2016-2017 годах вашему покорному слуге удалось съездить на полугодовое обучение в далекий город Ильменау (Германия), где он успешно (в общем и целом) закончил магистерскую программу Communications and Signal processing. Начну своё вступление издалека. Иногда... Программа оказалась не из простых, однако сейчас о ней вспоминать даже приятно.

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

Один из таких материалов перед вами.

Какие цели я преследовал, пока готовил семинар:

  1. рассказать о некоторых уже устоявшихся, "умных" подходах в теме антенных решеток наиболее доступно и сделать это на русском языке;
  2. провести небольшое моделирование на языке Python 3, чтобы сагитировать собратьев радиоинженеров присмотреться к языкам программирования (если ещё не присмотрелись);
  3. дать ссылки на хорошую англоязычную литературу — без чтения зарубежных источников, сейчас, увы, никуда.

Что рассмотрим:

  • Метод MUSIC (MUltiple SIgnal Classification) — к этому, собственно, и отсылка на превью.

Gist). Пример по диаграммообразованию и методу MVDR можно найти по ссылке (если по дополнительному материалу будут вопросы или предложения, то обсуждение можно продолжить на Github.

Как я уже сказал выше, использовать мы будем Python, а именно:

import numpy as np
import matplotlib.pyplot as plt

Потому что, мне хочется показать, что подобную работу можно сделать и на Python, а сфера применения Python куда шире, чем у MATLAB. Почему не MATLAB, один из самых популярных и удобных кандидатов для моделирования с использованием линейной алгебры, спросите вы? Поэтому быть знакомым с синтаксисом Python — дело, на мой взгляд, полезное.

Начнем!

Спасибо создателям за прекрасный инструмент! Формулы подготовлены через https://upmath.me/.

Постановка задачи

Допустим, имеется линейная антенная решетка, состоящая из некоторого количества элементов, разнесенных друг от друга на расстояние \Delta = \frac{2} (шаг антенной решетки), где \lambda — длинна несущей электромагнитной (ЭМ) волны.

На данную антенную решетку с разных направлений падают электромагнитные волны.

1. Рис. Адаптивная антенная система.

Как видно из рисунка, антенная решетка рассматривается, как адаптивный фильтр.

Собственно говоря, нахождение оптимального вектора коэффициентов (\mathbf{w}_{opt}) и есть основная задача адаптивных антенных решеток с математической точки зрения.

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

Моделирование принятого сигнала

Модель принятого сигнала мы можем представить через формулу:

\mathbf{X} = \mathbf{A} \mathbf{S} + \mathbf{N}

\quad \mathbf{a}(\theta_d)]"/> — матрица сканирующих векторов (steering vectors) антенной решетки (a_i = \exp(-j \mu m_i), <img src="https://habrastorage.org/getpro/habr/post_images/299/9cb/d24/2999cbd24c6ac4ba1fd3a74fdd0f5635.svg" alt="m = 0, 1 ... где \mathbf{A} = [\mathbf{a}(\theta_1) \quad \mathbf{a}(\theta_2) \quad ...  (M-1),M — количество элементов антенной решетки,d — количество источников ЭМ волн, \theta — угол направление прихода ЭМ волны), \mathbf{S} — матрица передаваемых символов, а \mathbf{N} — матрица аддитивных шумов.

ULA

2. Рис. 32]. Ненаправленная линейная антенная решетка (ULAA — uniform linear anntenna array) [1, с.

Мы не получаем в явном виде информацию о количестве источников и направлениях, однако, информация об этом всё же содержится в принятом сигнале. Переосмыслим данную формулу в "бытовом" ключе: на нашу решетку мы получаем некоторую "кашу" из различных сигналов, которую мы обозначаем через \mathbf{X}.

Начинаем искать!

Для этого обычно переходят к манипуляциям не с самими матрицами комплексных амплитуд сигналов, а с их ковариациями (т.е., по сути, с мощностями):

\mathbf{R}_{xx} = \mathbf{X}\mathbf{X}^H = \mathbf{A} \mathbf{R}_{ss} \mathbf{A}^H + \mathbf{R}_{nn}

Условия

Введем важное условие к рассмотрению: ограничение углового распознавания по Рэлею (Rayleigh angle resolution limit):

sin (\theta_R) = \frac{\lambda}{D}

где D = M \Delta — это длинна линейной решетки.

Переопределим угол прихода электромагнитной волны через понятие пространственной частоты:

\mu_R = \frac{2 \pi}{\lambda}\Delta sin(\theta_R) = \frac{2 \pi}{\lambda}\Delta \frac{\lambda}{\Delta M} = \frac{2 \pi}{M}

где \mu_R — есть стандартная ширина главного лепестка ДН (standard beamwidth).

Чтобы проверить насколько наш метод эффективен и при каких условиях, введем некоторые заданные значения для углового разделения:

  1. \mu_1 = - \mu_R, \quad \mu_2 = 0, \quad \mu_3 = \mu_R \quad — разделение в одну ширину луча;

  2. 5 \mu_R, \quad \mu_2 = 0, \quad \mu_3 = 0. \mu_1 = -0. 5 \mu_R \quad — разделение в одну вторую ширины луча;

  3. 3 \mu_R, \quad \mu_2 = 0, \quad \mu_3 = 0. \mu_1 = -0. 3 \mu_R \quad — разделение в три десятых ширины луча.

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

M = 10 # количество элементов решетки (сенсоров)
SNR = 10 # Отношение сигнал-шум (dB) d = 3 # количество источников ЭМ волн
N = 50 # количество "замеров" (snapshots) S = ( np.sign(np.random.randn(d,N)) + 1j * np.sign(np.random.randn(d,N)) ) / np.sqrt(2) # QPSK W = ( np.random.randn(M,N) + 1j * np.random.randn(M,N) ) / np.sqrt(2) * 10**(-SNR/20) # AWGN
# Общая формула:
# sqrt(N0/2)*(G1 + jG2), # где G1 и G2 - независимые гауссовские процессы.
# т.к. Es(энергия символа) для QPSK равна 1 Вт, спектральная мощность шума (noise spectral density): # N0 = (Es/N)^(-1) = SNR^(-1) [Вт] (принимаем в данном примере, что SNR = Es/N0); # или в логарифмическом масштабе:
# SNR_dB = 10log10(SNR) => N0_dB = -10log10(SNR) = -SNR_dB [дБ]; # Нам дано значение SNR в логарифмической шкале (т.е. в дБ), переводим в линейную:
# SNR = 10^(SNR_dB/10) => sqrt(N0) = (10^(-SNR_dB/10))^(1/2) = 10^(-SNR_dB/20) mu_R = 2*np.pi / M

Немного теории о самом методе

Рассматриваемая задача метода Писаренко заключалась в оценке частот суммы комплексных экспонент в белом шуме. Прежде всего отметим, что прародителем метода MUSIC является метод Писаренко (1973). Ф. В. Впоследствии этот метод стал особым случаем метода MUSIC. Писаренко продемонстрировал, что частоты можно найти из собственных векторов, соответствующих минимальному собственному значению автокорреляционной матрицы. 459] [2, c.

Основным подходом этого алгоритма является разложение матрицы ковариации принятого сигнала на собственные значения. Шмидт со своими коллегами предложил алгоритм классификации множественных сигналов (MUSIC) в 1979 году [4]. Здесь подпространства сигнала и шума вычисляются с использованием линейной алгебры, и являются ортогональными друг другу. Поскольку этот алгоритм учитывает некоррелированный шум, порожденная ковариационная матрица имеет диагональный вид. Поэтому алгоритм использует свойство ортогональности для выделения сигнальных и шумовых подпространств [5].

Обобщенный алгоритм MUSIC можно определить следующим образом:

  • Найти ковариационную матрицу \mathbf{R}_{xx}
  • Найти собственные вектора через EVD или же другой подходящий числовой алгоритм:

\mathbf{R}_{xx} = \mathbf{U}\mathbf{\Lambda}\mathbf{U}^H \qquad(1)

  • Найти псевдоспектр (почему с приставкой псевдо, мы обсудим ниже) MUSIC через следующую формулу:

P_{MU}(e^{j\omega}) = \frac{1}{\sum \limits_{i=d+1}^{M}|\mathbf{a}^H\mathbf{u}_i|^2} \qquad(2)

& e^{j(M-1)\omega} \end{bmatrix}^T"/> — вектор экспонент для частоты ω, лежащей в некотором заданном диапазоне, а \mathbf{u}_i — i-ый собственный вектор (eigen vector) ковариационной матрицы (1), соответствующие шумовому подпространству матрицы (1) — отсюда и индексация с d+1 (d — ранг матрицы (1)). где <img src="https://habrastorage.org/getpro/habr/post_images/bb7/7ab/1aa/bb77ab1aa0a09b1807e947dcd80b40c8.svg" alt="\mathbf{a} = \begin{bmatrix} e^{j0\omega} & e^{j1\omega} & e^{j2\omega} & ...

Обратите внимание на два основных момента:
Для большей наглядности попробуйте прогнать соответствующий MATLAB-скрипт представленный по ссылке.

  • вместо вычисления квадрата второй нормы в знаменателе (2) авторы применяют к собственным векторам алгоритм БПФ, что облегчает моделирование за счет применения встроенных функций и, в общем, не противоречит теории с математической точки зрения;
  • вычисление ковариационной матрицы производится через сверточные матрицы, выше для оценки пространственных частот был показан другой подход.

Алгоритм вычисления псевдоспектров в данном контексте приведем ниже: Как нетрудно догадаться из названия, MUSIC также является классическим методом оценки направления приема с высоким разрешением.

\mathbf{U} = [\mathbf{U}_s \quad \mathbf{U}_0]

  • выбираем некоторый диапазон поиска:

a(\mu) = \begin{bmatrix} e^{j0\mu_1} & ... & e^{j0\mu_Q} \\ ... &...&... \\ e^{j(M-1)\mu_1} & ... & e^{j(M-1)\mu_Q} \end{bmatrix}

где \mu = -\frac{2\pi f_c}{c}\Delta sin\theta = -\frac{2\pi }{\lambda }\Delta sin\theta

  • вычисляем псевдоспектр:

P_{MU}(\theta)=\frac{\mathbf{a}^H (\theta)\mathbf{a}(\theta)}{\mathbf{a}^H(\theta)\mathbf{U}_0 \mathbf{U}_0^H \mathbf{a}(\theta)}

Связь между спектральным анализом и анализом углов прихода (DoA — direction of arriaval) ЭМ волн описана в таблице 1.

Таблица 1 Связь между приложениями MUSIC: Обработка массива сигналов и Гармонический поиск [6].

Переменная

Обработка массива сигналов

Гармонический поиск

M

Количество сенсоров

Количество временных отрезков

N

Количество временных отрезков

Количество экспериментов

d

Количество волновых фронтов

Количество комплексных компонент

\mu

Пространственные частоты

Нормализованные частоты

в сущности каждый сенсор, принимая волну с определенной задержкой фазы (т.е. В общем и целом, процесс приема через массивы (решетки) можно сравнить с процессом классической дискретизации, т.к. Количеству реализаций (экспериментов) классического спектрального анализа будет соответствовать количество временных отрезков (snapshots). с определенной временной задержкой), выполняет функции дискретизирующего дельта импульса. Каждому источнику будет соответствовать свой волновой фронт, что эквивалентно количеству уникальных синусоид сигнала в случае спектрального анализа.

Мы уже упомянули выше, что вектора a(\theta_i)\epsilon A, где i=1,2,..,d ортогональны шумовому подпространству ковариационной матрицы, т.е.: А теперь вернемся к моменту вычисления собственных векторов.

a(\theta_i)^TU_0=0^T

Такой метод в отличие от числовых алгоритмов (к которым, как мы отметили выше, относится и EVD) позволяет получить настоящие, а не приближенные собственные значения. Собственно говоря, мы видим систему уравнений, решая которую мы можем найти корни — собственные вектора. Эта же идея легла в основу алгоритма Root MUSIC. Именно поэтому данный подход позволяет получить не псевдоспектр, а спектр.

Моделирование

Наконец-то все формулы описаны и сколько-то объяснены. Фуф! Можем приступать к моделированию.

cases = [[-1., 0, 1.], [-0.5, 0, 0.5], [-0.3, 0, 0.3],]
for idxm, c in enumerate(cases): # углы прихода (пространственные частоты): mu_1 = c[0]*mu_R mu_2 = c[1]*mu_R mu_3 = c[2]*mu_R # сканирующие вектора a_1 = np.exp(1j*mu_1*np.arange(M)) a_2 = np.exp(1j*mu_2*np.arange(M)) a_3 = np.exp(1j*mu_3*np.arange(M)) A = (np.array([a_1, a_2, a_3])).T # матрица сканирующих векторов X = np.dot(A,S) + W # матрица принятых сигналов R = np.dot(X,np.matrix(X).H) U, Sigma, Vh = np.linalg.svd(X, full_matrices=True) U_0 = U[:,d:] # шумовое подпространство thetas = np.arange(-90,91)*(np.pi/180) # диапазон азимутов mus = np.pi*np.sin(thetas) # диапазон пространственных частот a = np.empty((M, len(thetas)), dtype = complex) for idx, mu in enumerate(mus): a[:,idx] = np.exp(1j*mu*np.arange(M)) # MVDR: S_MVDR = np.empty(len(thetas), dtype = complex) for idx in range(np.shape(a)[1]): a_idx = (a[:, idx]).reshape((M, 1)) S_MVDR[idx] = 1 / (np.dot(np.matrix(a_idx).H, np.dot(np.linalg.pinv(R),a_idx))) # MUSIC: S_MUSIC = np.empty(len(thetas), dtype = complex) for idx in range(np.shape(a)[1]): a_idx = (a[:, idx]).reshape((M, 1)) S_MUSIC[idx] = np.dot(np.matrix(a_idx).H,a_idx)\ / (np.dot(np.matrix(a_idx).H, np.dot(U_0,np.dot(np.matrix(U_0).H,a_idx)))) plt.subplots(figsize=(10, 5), dpi=150) plt.semilogy(thetas*(180/np.pi), np.real( (S_MVDR / max(S_MVDR))), color='green', label='MVDR') plt.semilogy(thetas*(180/np.pi), np.real((S_MUSIC/ max(S_MUSIC))), color='red', label='MUSIC') plt.grid(color='r', linestyle='-', linewidth=0.2) plt.xlabel('Azimuth angles θ (degrees)') plt.ylabel('Power (pseudo)spectrum (normalized)') plt.legend() plt.title('Case #'+str(idxm+1)) plt.show()



Как мы можем заметить, MUSIC обладает большей разрешающей способностью и позволяет добиться, в целом, лучших результатов, чем, например, позволяет MVDR — такой же представитель параметрических методов спектрального анализа.

Однако нужно учитывать, что при использовании MUSIC мы задействуем более дорогие в вычислительном плане алгоритмы, такие как EVD или SVD, что является некоторой ценой за большую точность.

Такие дела.

Список использованной литературы:

  1. Haykin, Simon, and KJ Ray Liu. Handbook on array processing and sensor networks. Vol. 63. John Wiley & Sons, 2010. pp. 102-107
  2. Hayes M. H. Statistical digital signal processing and modeling. – John Wiley & Sons, 2009.
  3. Haykin, Simon S. Adaptive filter theory. Pearson Education India, 2008. pp. 422-427
  4. Richmond, Christ D. "Capon algorithm mean-squared error threshold SNR prediction and probability of resolution." IEEE Transactions on Signal Processing 53.8 (2005): 2748-2764.
  5. S. K. P. Gupta, MUSIC and improved MUSIC algorithm to esimate dorection of arrival, IEEE, 2015.
  6. Professor Martin Haardt's lectures (array array)
Показать больше

Похожие публикации

Добавить комментарий

Ваш адрес email не будет опубликован. Обязательные поля помечены *

Кнопка «Наверх»