Алгоритм MUSIC: MUltiple SIgnal Classification в оценке спектра

Вот мы наконец и добрались до этого метода. Литературы по этой технике — целая куча, но как обычно смысл ускользает. Само описание MUSIC достаточно простое: это поиск экстремума произведения двух матриц. Однако что происходит внутри этих матриц — с этим нам предстоит разобраться.

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

Будем рассматривать один частный случай — пример входного сигнала для алгоритма, и на основе этого случая будем делать обобщения на все возможные случаи. За доказательством справедливости этих обобщений отсылаем дотошных читателей к многочисленной литературе; здесь во главу угла поставлена простота и доходчивость изложения.

Компактный симпатичный спектр из двух частот

Зададим исходные данные. Будем иметь дело с сигналом, состоящим из двух частот; также подмешаем немного шума, для реальности происходящего:

На самом деле, мы будем наблюдать смесь этих сигналов, и наша задача — определить какие частоты присутствуют в его спектре. Поэтому наблюдаемый сигнал X будет суммой этих массивов:

Суммирование идет по выборкам, о чем говорит выбор оси axis=1. Поэтому размерность полученной матрицы будет N x N.

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

А теперь — прощай, сигнал X! В нашем алгоритме ты нам больше не понадобишься. Вместо тебя будет работать твоя корреляционная матрица Rxx. Сразу скажу, что быть на сцене ей тоже осталось недолго: после того как мы находим собственные числа en и векторы ev корреляционной матрицы Rxx, она тоже перестает быть нужной:

Гарантирую, что с этого момента вы потеряли нить и вам расхотелось читать дальше ) Поэтому самое время остановиться и понять смысл того, что мы делаем. Тем более что большую часть алгоритма мы уже прошли: осталось фактически только одно последнее действие ) Вот она, обманчивая простота матриц!

Итак, наступил обещанный момент погружения в матричные пространства.

Зачем нужны собственные числа и векторы

Идея метода MUSIC состоит в разделении пространства входных сигналов на сигнальную часть и помеховую. Для пространства сигналов существует свой базис, в котором может быть разложена корреляционная матрица Rxx. Поскольку Rxx является самосопряженной (сравниваем Rxx и Rxx.H в Питоне — это одно и то же), то базис разложения будет ортогональным и будет представлять из себя не что иное, как собственные векторы матрицы, масштабированные собственными числами. Что такое ортогональное разложение вы конечно же помните из школьной геометрии, когда рисовали проекции вектора на ортогональные оси x и y. Только это было двумерное пространство, которое легко представить геометрически, а в нашем случае размерности пространства N только и остается, что положиться на формулы.

Проверим, соответствует ли Rxx разложение по ортогональному базису:

результат будет соответствовать Rxx.

Теперь мы можем предположить, что одна часть собственных векторов, участвующих в разложении Rxx, соответствует сигналу, а другая часть — помехе. Будем также держаться предположения, что полезных сигналов — P, а все остальное в количестве N-P это помехи или шум. Наше предположение некоторым образом подтвердят собственные числа, если мы расставим их в порядке убывания:

По собственным числам можно судить об уровне входных сигналов: как и ожидалось, первые P=2 это полезный гармонический сигнал, остальное — это помеха. Не откладывая на потом, сразу упорядочим и собственные векторы, чтобы они соответствовали собственным числам, расположенным по убыванию:

Теперь выстраиваем следующую логическую цепочку. Начало ей положит тот факт, что собственные векторы сигнала и помехи — взаимно ортогональны:

результат — вектор практически близкий к нулю, что и означает ортогональность.

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

Разложение экспоненциального вектора на базис помехи

Теперь все дело за тем, как правильно записать нотацию искомого разложения. Собственные векторы помехи мы уже выделили раньше: это вектор EV. Вектор комплексных опорных сигналов будет иметь размерность N:

Значения частот будем перебирать из диапазона:

Абсолютное значение спектра будем находить как модуль разложения REF*EV.H*EV*REF.H.

Тогда полный цикл перебора опорного вектора будет выглядеть так:

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

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

показан на рисунке.

Начнем с P=2:

MUltiple SIgnal Classification, MUSIC, Eigenvectors, Eigenvalues, Оценка спектра, Корреляционная матрица, Correlation Matrix

Оценка спектра MUSIC для P=2

Здесь неожиданностей нет. На рисунке четко видно два пика относительной частоты 3.0 и 5.0, причем уровень второго сигнала ниже. Мы четко определили наличие двух синусоид в комплексном сигнале в присутствии шума, причем получили значение частот этих сигналов и относительные амплитуды. Теперь возникает интересный вопрос: что произойдет если мы зададим P=1? По идее алгоритм должен посчитать за полезный только один сигнал, а вторую синусоиду вместе с шумом отнести к помехе. Так оно и есть:

MUltiple SIgnal Classification, MUSIC, Eigenvectors, Eigenvalues, Оценка спектра, Корреляционная матрица, Correlation Matrix

Оценка спектра MUSIC для P=1

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

Продолжаем эксперимент.  Теперь предположим, что на входе три полезных сигнала, то есть P=3. Что произойдет если белый шум также будет засчитан на равных как нечто полезное?

MUltiple SIgnal Classification, MUSIC, Eigenvectors, Eigenvalues, Оценка спектра, Корреляционная матрица, Correlation Matrix

Оценка спектра MUSIC для P=3

Получилось довольно неприятно, но не смертельно. Третий всплеск конечно достаточно большой по амплитуде, и такой же всплеск наблюдаем на нулевой частоте. Но мы сами задали шуму приличный уровень — 0.5, так что метод MUSIC довольно лояльно относится к тому, что мы неточно определили количество полезных сигналов.

Продолжим расширять наше представление о количестве полезных компонент и зададим P=4:

MUltiple SIgnal Classification, MUSIC, Eigenvectors, Eigenvalues, Оценка спектра, Корреляционная матрица, Correlation Matrix

Оценка спектра MUSIC для P=4

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

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

 

 

 

Ответить

Вы можете использовать эти HTML теги

<a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code class="" title="" data-url=""> <del datetime=""> <em> <i> <q cite=""> <s> <strike> <strong> <pre class="" title="" data-url=""> <span class="" title="" data-url="">