CC BY
97
2007
НАУЧНЫЙ ВЕСТНИК МГТУ ГА серия Радиофизика и радиотехника
№126
УДК 621.396
НЕКОТОРЫЕ АЛГОРИТМЫ РАЗВЕРТЫВАНИЯ ФАЗОВОГО СПЕКТРА
Ю.А. КРАСНИТСКИЙ
Рассматриваются несколько широко распространенных алгоритмов вычисления развернутого фазового спектра. В качестве меры несовпадения результатов их применения к одному и тому же одномерному сигналу предложено использовать величину угла между векторами этих спектров, вычисляемую через скалярное произведение последних. Роль опорного отведена алгоритму, который работает на основе анализа расположения нулей 2-преобразования исследуемого сигнала вблизи единичной окружности и обладает наибольшей устойчивостью. Показано, что основная причина расхождений состоит в выборе недостаточно малого шага взятия отсчетов в частотной области.
1. Введение
При исследовании сигналов различной физической природы часто возникает задача нахождения их фазовых спектров. Она, в частности, представляет собой важный этап вычисления комплексного кепстра, а также нахождения зависимости мгновенной частоты изучаемого сигнала от времени на основе преобразования Г ильберта.
Фазовый спектр сигнала и(/) может быть получен следующим образом.
1) Найдем комплексный спектр
£(т) ехр [і аг§ £(т)]
(1)
£ (/о) = Б1 \и(г)]
этого сигнала. Символом Б1 обозначено прямое преобразование Фурье. Стандартный алгоритм вычисления аргумента указанного преобразования имеет вид
аг§ £ (ш) = аг^ [1т £ (Ш)/Яе £ (/о)] (2)
и дает главные значения фазы сигнала на частотах со, определенные по модулю 2р. Они лежат в промежутке \-р р], на границах которого функция (2) в случае, если сигнал и{() не является минимально фазовым, претерпевает разрывы.
Это подтверждает рис. 1, где показаны два сигнала различной физической природы - атмо-
сферик (рис. 1 а), порожденный излучением молниевого разряда, и сигнал на выходе приемника георадара (рис. 1с). Фазовые спектры (Ь и ф вычислены с помощью (2). Разрывы затрудняют интерпретацию этих спектров и делают их неудобными для дальнейшего использования.
2000
20
60
100 120 0 Рис. 1
200 400 600 800 1000
2) С другой стороны, комплексный логарифм спектра (1), равный
(іт) = Ьо§ £ (т) = 1п
£(ш) + /[аг^(ш) + 2ят^ , т = 0,± 1,±2... (3)
есть бесконечнозначная функция. Ее действительная часть определяется однозначно, а мнимая
<р(о) = 1т (/о) = ащ £ (/о) + 2рт , т = 0, ± 1, ± 2 . (4)
которая и есть фазовый спектр, - с точностью до слагаемого, кратного 2р.
Поскольку комплексный логарифм (3) должен быть непрерывной функцией, его мнимую часть, т.е. фазовый спектр (4), нельзя ограничивать пределами \-р р], в которых лежит главное
значение последнего, соответствующее т = 0. Операция устранения многозначности в (4) и (3) состоит в восстановлении непрерывной зависимости рассматриваемых функций от частоты и называется развертыванием фазы.
При цифровой обработке одномерных сигналов применяют около 10 различных алгоритмов развертывания фазы и их модификаций. Они часто основаны на эвристических соображениях, сложны в использовании и не обладают устойчивостью по отношению к небольшим вариациям главного значения фазы и малым изменениям шага дискретизации по частоте. Рассмотрим наиболее употребительные алгоритмы, предполагая, что анализируемый сигнал задан отсчетами во временной области, а спектры (1) - (4) представляют собой дискретные функции частоты.
2.1. Алгоритм Шафера. Он базируется на том, что к выборкам главного значения фазы, которые найдены в (2), добавляют величины, равные 2рт, где т - целое число. Значения т образуют так называемую корректирующую последовательность и определяются путем сравнения соседних отсчетов фазы. Метод подробно описан в [1,2] и реализован в большинстве пакетов программ обработки сигналов, в том числе и в МайаЬ [3]. Главный его дефект состоит в том, что он неустойчив (в смысле, оговоренном выше) при недостаточно малом шаге взятия частотных выборок в (2).
2.2. Алгоритм Триболе основан на адаптивном интегрировании производной фазового спектра по частоте и может быть реализован в результате следующей последовательности действий.
а) Находят спектр
обрабатываемого дискретного сигнала и(п), записав (2) в алгебраической форме, находят спектр
дискретного сигнала пи(п). Для получения спектров (5) и (6) обычно применяют быстрое преобразование Фурье.
б) Вычисляют производную фазового спектра, используя соотношение
в) Получают развернутый фазовый спектр путем адаптивного интегрирования производной (7), вычисляя при необходимости ее дополнительные отсчеты, т.е. интерполируя в промежутке \п п+1] с помощью дискретного преобразования Фурье.
Рассмотренный алгоритм подробно описан в [4]. Им пользуются при обработке сейсмических сигналов. Существенный его недостаток заключается в росте времени, затрачиваемого на обработку, за счет применения интерполяции производной (7).
2.3. Применение медианной фильтрации. Для решения задачи развертывания фазового спектра может быть использовано известное свойство медианного фильтра. Если на его вход подать соответствующую последовательность отсчетов, в состав которой входят прямоугольные импульсы шириной к (т.е группы, содержащие по к отсчетов одинаковой амплитуды), то на выходе фильтра эти импульсы будут подавлены. Порядок фильтра составляет 2к + 1, где к - целое положительное число. Таким образом, для подавления одиночных импульсов шириной в один отсчет требуется медианный фильтр 3-го порядка. Если между двумя уровнями сигнала существует линейно нарастающая (или спадающая) переходная область, то для ее устранения путем медианной фильтрации сигнал следует предварительно продифференцировать, превратив область нарастания или спада в прямоугольный импульс.
2. Некоторые методы вычисления развернутого фазового спектра
$,(п) = Б1 [и(п)] = Яе $>(п) + і 1т $,(п)
(5)
81(п) = Б1 [пи(п)] = Яе 81(п)+і 1т 81(п) .
(6)
ф'п) = Яе 5(п)Яе $(п) + 1т $>(п) 1т $(п)
(7)
На основе этих соображений разработан алгоритм развертывания фазы, описанный в [5].
Он учитывает, что разрывы непрерывности в неразвернутом фазовом спектре (3) имеют еди-
ничную ширину, и включает в себя следующие операции.
a) Вычисление разностей
Ф'(n) = arg S(n) - arg S(n -1) (В)
соседних отсчетов неразвернутого фазового спектра (4).
b) Нахождение корректирующей добавки
АФ(п) = median { [Ф'(n -1) Ф'(n) Ф'(n +1)]} (9)
к фазе путем обработки вектора, составленного из трех последовательных значений функции (В), фильтром 3-го порядка, который оценивает скользящую медиану.
c) Рекурсивное вычисление развернутой фазы
Ф0(п) = arg S (n), n = 1,
Ф0(п) = Ф0(п -1) + АФ(п), n = 2,3...
Следует отметить, что сведения о подробных исследованиях рассмотренного способа получения непрерывного фазового спектра в доступных источниках отсутствуют.
2.4. Рекурсивный алгоритм Крайника отличается простым и изящным способом оценки приращения ЛФ(п) в (10), основанным [б] на вычислении аргумента отношения соседних отсчетов спектра (1). Действительно,
(10)
S(n +1) | S(n +1) | єхр[/Ф(п +1)]
arg —----- = argJ!----------———— -----—
= Ф(п +1) - Ф(п) = А Ф(п) .
(11)
8(п) 18(п) | ехр[/Ф(п)]
В отличие от алгоритма Шафера здесь для нахождения развернутого фазового спектра не требуется ни обнаружения разрывов, ни добавления величин, кратных 2р, к главным значениям неразвернутой фазы (2). Нет также необходимости проводить интерполяцию между соседними отсчетами, как в алгоритме Триболе.
2.5. Развертывание фазы на основе факторизации спектра было предложено в 1982 г. [7], но из-за больших вычислительных затрат долгое время не применялось на практике. В этом методе используется анализ расположения нулей 2-преобразования исследуемого сигнала, лежащих в узкой полосе вблизи единичной окружности. На рис. 2 в качестве примера показаны нули 2-преобразования атмосферика, изображенного на рис. 1а.
1.5
1
0.5
0
•
•
" • (.
-10 1;
Рис. 2
Если нуль кратности 1 располагается произвольно близко к единичной окружности и находится между частотами сои со+Лсо, то скачок главного значения фазы в (2) составит ±р Выбор знака зависит от того, лежит нуль вне (+) или внутри (-) этой окружности. При расположении нуля непосредственно на окружности величина скачка будет равна 2р
Практическая реализация этого метода развертывания фазы стала возможной лишь после разработки эффективного алгоритма факторизации полиномов высоких степеней [8].
3. Сравнение методов
В [8] указывается, что применение алгоритма, основанного на факторизации спектра, гарантирует точное решение задачи развертывания фазы. Этот результат может быть принят за эталон при тестировании других методов нахождения развернутого фазового спектра.
В качестве меры близости к эталону выберем величину угла поворота вектора развернутой фазы, вычисленного по одному из описанных выше алгоритмов 2.1, 2.3 или 2.4, относительно того же вектора, найденного по методу факторизации 2.5. Этот угол можно определить из соотношения
соб А =
Ф, Ф,
уо
Ф
Фо
уо
(12)
где || Ф|| и || Фэт|| - энергетические нормы, а угловые скобки обозначают скалярное произведение
(13)
сравниваемого и эталонного векторов. Здесь Юм - частота Найквиста.
Цель сравнения алгоритмов состоит в том, чтобы выяснить, как зависит ошибка в развертывании фазового спектра от величины шага в частотной области и как этот шаг влияет на время вычислений.
На рис. 3 представлены результаты численного эксперимента, состоявшего в том, что при вычислении фазового спектра атмосферика, представленного на рис. 1 а, шаг по частоте мог быть последовательно уменьшен путем добавления к исходному сигналу нулевых отсчетов. Их число выбиралось таким, чтобы длина обрабатываемого сигнала стала равной 2п, где п = 8 ...
11. Такой же размер имел формат быстрого преобразования Фурье, применявшегося при вычислениях. На рис. 3 он обозначен через КУ.
0
Рис. 3
Графики неразвернутых и развернутых фазовых спектров показаны при изменении частоты от 0 до частоты Найквиста. Кривая 1 обозначает фазовый спектр, полученный на основе алгоритма факторизации и принятый в (12) и (13) за эталон. Цифрами 2, 3 и 4 помечены результаты развертывания фазы по алгоритмам Шафера, Крайника и с использованием медианной фильтрации соответственно.
Анализ приведенных кривых показывает, что по мере уменьшения шага по частоте, что достигается путем увеличения размерности КУ преобразования Фурье, они сближаются. Для получения приемлемых результатов число дополнительных нулевых отсчетов должно быть таким, чтобы исходная длина сигнала увеличилась бы не менее, чем в 4 раза. Наибольшее расхождение относительно эталона характерно для спектра, вычисленного по методу медианной фильтрации.
К таким же выводам приводит и анализ поведения углов (12) поворота векторов развернутой фазы. На рис. 4 приведены результаты их расчета (в градусах) для обоих сигналов,
показанных на рис. 1. Через п обозначены показатели степени, в которые надо возвести 2, чтобы найти длину сигнала с учетом добавления нулей.
Таблица1
Быстродействие (сек) алгоритмов развертывания фазы
n 1 2 3 4
11 7.53 0.05 0.66 2.91
12 7.41 0.11 1.97 6.32
13 7.19 0.11 6.48 17.03
14 6.86 0.17 20.37 37.6
15 6.69 0.27 73.1 107.8
Рис. 4
Для оценки быстродействия сравниваемых алгоритмов была применена функция cputime системы Matlab, которая позволяет найти вычислить процессорное время. Результаты этого исследования на примере сигнала s070ac (рис. 1с) показаны в табл. 1, где колонки с 1 по 4 обозначены так же, как и кривые на рис. 3. Вычисления выполнены на компьютере с процессором Pentium (r) II Intel MMX (TM).
4. Выводы
Проведенные исследования подтверждают, что характеристики алгоритмов развертывания фазы сильно зависят от шага в частотной области. Приемлемые для большинства приложений результаты достигаются при не менее чем 4-кратном удлинении обрабатываемого сигнала путем добавления нулевых отсчетов. Точностные характеристики алгоритмов Шафера и Крайника практически одинаковы, но быстродействие первого из них существенно выше. Наихудшие показатели по точности и быстродействию демонстрирует алгоритм, использующий медианную фильтрацию.
ЛИТЕРАТУРА
1. Oppenheim A.V., Schafer R.W. Discrete-Time Signal Processing. - N.-Y.: Englewood Cliffs, Prentice Hall,
1989.
2. Васильев В., Гуров И. Компьютерная обработка сигналов в приложении к интерферометрическим системам. - СПб.: БХВ-Санкт-Петербург, 1998.
3. Рудаков П.И., Сафонов В.И. Обработка сигналов и изображений. - М.: ДИАЛОГ-МИФИ, 2000.
4. Tribolet J. M. Seismic Applications of Homomorphic Signal Processing. - N.-Y.: Englewood Cliffs, Prentice Hall,
1977.
5. Loeffler C.M., R.E. Leonard, Jr. Phase unwrapping via median filtering. Int. Conf. on Acoust., Speech and Sign. Proc. (ICASSP84), 1984, 46.8/1 - 46.8/3 - 1982,
6. Krajnik Е. A simple and reliable phase unwrapping algorithm. Signal Processing VI: Theories and Applications. J. Vandewalle (eds.). - Elsevier, 1992, 917-919.
7. Steiglitz K., Dickinson B. Phase unwrapping by factorization. IEEE Trans. on Acoust., Speech and Sign. Proc., 1982, v.30, No 6, 984-991.
8. Fox J.W., Sitton G.A., Lindsey J.P., Burrus C.S., Treitel S. Factoring ultra - high degree polinomials. - Rice University, Nov. 26, 2003 (www-dsp.rice.edu/software).
SOME ALGORITHMS OF PHASE SPECTRUM DEPLOYMENT
Krasnitsky Yu.A.
Some widespread algorithms of calculation of deployed phase spectrum as considered. As mismatch measure of their application results to the same regular signal it is proposed to use the angle value between vectors of these spectra which is defined by one is scalar product. The reference algorithm is that based on the analysis of zero pole location of z-transformation of analyzed signal near to unit circumference and it has maximum stability. It is shoran that the main reason of variants is in selection of not enough short stop of sampling in frequency domain.
Сведения об авторе
Краснитский Юрий Александрович, 1937 г.р., окончил Ленинградский институт точной механики и оптики (1960), доктор технических наук, профессор РАИ, автор более 100 научных работ, область научных интересов - антенные устройства, обработка сигналов, однопунктная радионавигация.
