CC BY
108
УДК 004.67
Н.Т. Сафиуллин, С.В. Поршнев
сравнительный анализ расчета мгновенной частоты через преобразование гильберта и прямую квадратуру
Понятие аналитического сигнала (АС), предложенное в 1946 г. Габором [1], позволило ввести определения амплитуды, фазы и мгновенной частоты сигнала м(?). По Габору АС м>(£), соответствующий действительной функции м(?), определяется для всех функций, принадлежащих классу
L р (-да, ,
м>Ц) = м(?) + рЦ) = а(Г)вм,), (1)
где у({) — сигнал, сопряженный к исходному сигналу по Гильберту,
v(t) =1 +Г ^ dт = H (u ), п „ t - т
(2)
здесь H(u) - оператор, описывающий преобразование Гильберта; a(t) — функция, имеющая физический смысл огибающей сигнала u(t),
a(t) = |w(t) = ylu2(t) + v2(t) . (3)
Использование АС, определенного в соответствии с (1), позволяет естественным путем ввести для сигнала u(t) понятие фазы ip(t)
9(t) = arctg I v(t) | = и (t )
= arccos
u(t) 1 • f v(t) —^ | = arcsin | a(t ) J ^ a(t )
и мгновенной частоты ra(t) сигнала
u(t )v(t ) - u(t)v(t )
ra(t) = ф (t) = -
a1 (t)
(4)
(5)
После опубликования работы [1] на протяжении нескольких десятилетий вопросы, связанные с различными аспектами применения ПГ и трактовки характеристик сигналов, вычисляемых в соответствии с (1)—(5), стали предметом многочисленных научных исследований. При этом их авторы высказывали диаметрально противоположные мнения о физической сущности огибающей и мгновенной частоты сигнала, вводимых через ПГ, и правомерности использования данных величин для обработки сигналов и объяснения работы радиотехнических устройств [2—12]. Например, в [2] говорится, что «в общем случае Л(() и
©(О (огибающая и мгновенная фаза - прим. авт.), формально полученные посредством описанной процедуры, не будут медленными функциями и, называя их мгновенной амплитудой и фазой, мы не извлечем из этого никаких оправданных последствий». А в [3], в свою очередь, утверждается: «.. .любое определение и применение колебательных характеристик правильно, лишь постольку это определение и применение согласуется с тем, что дает аналитический сигнал. Иначе говоря, нет колебательных задач, в которых амплитуду, фазу и частоту следует вводить иначе, без аналитического сигнала». Отметим, что после выхода в свет монографии [3] произошло заметное снижение интереса исследователей к обсуждаемым проблемам, продолжавшееся вплоть до середины 90-х гг. ХХ в.
Возобновление исследований возможности использования ПГ в задачах обработки реальных сигналов было обусловлено необходимостью обработки сигналов и({), получаемых с выхода до-плеровского радиолокатора непрерывного излучения диапазона СВЧ при измерении параметров движения снаряда в стволах артиллерийских систем и стрелковых вооружений во время выстрела (типичные частоты электромагнитных волн (ЭВ), используемые в данных радиолокаторах, 10 ГГц и 33 ГГц) [4]:
4П - 1 (6)
и ^) = а({ ^т х^)
где X — длина ЭВ; а({) — функция, описывающая амплитудную модуляцию, обусловленную откатом ствола и воздействием случайных факторов; х(^ — перемещение боеприпаса в стволе.
Данные сигналы, зарегистрированные в цифровой форме на конечном временном интервале, представляют собой частотно-модулированные (ЧМ) сигналы, у которых зависимость мгновенной частоты от времени описывается некоторой монотонно возрастающей функцией. В зависимости от типа артиллерийской системы, а также от частот используемых ЭВ, длительности сигналов, полу-
чаемых в рассматриваемой задаче, варьируются в диапазоне от 1 до 10 мс, а полосы частот, занимаемые данными сигналами в спектральной области, - соответственно, в диапазонах [0,50] КГц до [0,400] КГц. Здесь обработка сигнала состоит в нахождении по значениям радиолокационного сигнала щ , измеренным в узлах дискретной временной сетки конечной длительности = 1, Ы, мгновенных значений фазы ) и частоты / ) сигнала и далее соответствующих перемещения и скорости движения боеприпаса:
х(*1) = ^Ф(<1),
(7)
4п
^) = ^ / й). 4п
В связи с тем, что рассматриваемые сигналы относятся к классу цифровых сигналов конечной длительности, т. е. не принадлежат пространству L р (-го, , потребовалось проведение целенаправленных исследований возможности использования ПГ для анализа данного класса сигналов. Полученные результаты, изложенные в [5], помогли сделать следующие выводы,
а также предложить соответствующие алгоритмы, позволяющие по функции ювыч(0 находить функцию ):
1) для широкополосных сигналов конечной длительности а(? )cos(ф(t)) понятие огибающей, определяемой через амплитуду соответствующего АС, является физически содержательным для сигналов как с узкополосной, так и с широкополосной несущей;
2) при перекрытии спектров функции амплитудной модуляции а(?) и несущей соб(ф(^)) вычисляемая огибающая сигнала авыч(?) имеет систематическую погрешность Аа(?) (авыч () = а(?) + Аа(?)), которая может быть уменьшена полосовой фильтрацией функции авыч(?) ;
3) причина искажения вычисляемой огибающей сигнала конечной длительности а)соб(ф(?)) - перекрытие спектров функций, описывающих закон амплитудной модуляции а(?) и несущей cos(ф(t));
4) в связи с тем, что аввыч(?) Ф а(?), вычисляемый закон изменения частоты широкополосного сигнала конечной длительности во времени
Рис. 1. Схема алгоритма, реализующего метод EMD
4
Рис. 2. Схема алгоритма, реализующего метод ПК (DQ)
ювыч(0 также имеет систематическую погрешность (ювыч (0 = ®(t) + Дш(?)).
ПГ вновь оказалось востребованным в 1998 г., когда Н. Гуанг предложил новый метод анализа нестационарных временных рядов, получивший название «преобразование Гуанга-Гильберта» (HHT) [6]. Данный метод основан на разделении исходного сигнала на квазигармонические компоненты c(t) с помощью алгоритма эмпирической модовой декомпозиции (EMD), блок-схема которого представлена на рис. 1. Далее, в методе ННТ производится расчет значений огибающей и мгновенной частоты каждой из выделенных компонент (или мод, как их чаще называют в иностранной литературе) в соответствии с (1)-(5).
Анализ опыта использования метода HHT для анализа реальных сигналов, имеющих конечную длительность, привел к тому, что результаты, имеющие однозначную трактовку с физической точки зрения, удается получить только на первом этапе реализации данного метода (т. е. при декомпозиции). В то время как при формальном применении ПГ к выделенным модам в ряде случаев значения МЧ противоречили физическим представлениям об изучаемых процессах [7].
В 2009 г. Н. Гуанг предложил новый способ расчета МЧ сигнала, в котором используется прямая квадратура (ПК или DQ - Direct Quadrature) [7]. В основе метода ПК, блок-схема которого представлена на рис. 2, лежит разложение исходного сигнала u(t) на амплитудную a(t) и фазовую ф(0 составляющие эмпирическим способом с дальнейшим расчетом МЧ в (5). Этот метод, по мнению его автора, позволяет получать результаты, намного превосходящие ПГ по точности. В связи с этим
представляет практический интерес проведение сравнения точностных характеристик цифровых ЧМ сигналов конечной длительности, получаемых с помощью ПГ и метода ПК.
Анализ результатов обработки ЧМ сигнала
Рассмотрим результаты обработки модели радиолокационного сигнала, частотно-временные характеристики которого подобны сигналу, получаемому при измерении параметров движения реального боеприпаса в стволе артиллерийской системы: длительность сигнала TS = 8,96 мс, полоса частот в спектральной области - [0; 30] КГц, частота дискретизации сигнала 7,3 МГц. Закон изменения частоты данного сигнала во времени представлен на рис. 3 а (график 1) (отметим, что здесь приходится использовать закон изменения частоты сигнала, заданный таблично в связи с отсутствием адекватно описывающих его аналитических выражений).
Результаты вычисления МЧ модели радиолокационного сигнала в соответствии с ПГ представлены на рис. 3 а (график 2). Отметим, что на практике вместо численного вычисления интеграла
T
= if sin(<P(т))
тт J
(8)
п Р г-т
оказывается удобным использовать известное соотношение между спектром исходного сигнала F(&) и спектром соответствующего ему АС Ж(ю) [3]:
2^(ш), ш > 0, Ж(ш) = ] ^(0), ш = 0, (9)
0, ш< 0.
a)
б)
х 10
, №1 __________
0 5 х 10 2 ПГ 6 8 * ю"3
0 х 10 2 ПК 6 8 х 10"
№3
•
х 10'
Рис. 3. Зависимости мгновенной частоты (а) и логарифма относительной погрешности вычисления МЧ (б)
модели радиолокационного сигнала от времени
График 3 демонстрирует результаты вычисления МЧ модели радиолокационного сигнала в соответствии с методом ПК.
Зависимость относительной погрешности вычисления МЧ модели радиолокационного сигнала каждым из методов
' ДО-со*» I
AcoJ0^
(10)
со(0
представлена на рис. 3 б, из которого видно, что для рассматриваемого сигнала, за исключением участков на краях интервала существования сигнала, где оба метода показали неудовлетворительную точность (относительная погрешность каждого может превышать 50 %), точность метода ПК оказалась на порядок выше ПГ. Результат, полученный с помощью ПГ, является вполне ожидаемым и объясняется наличием описанной выше систематической составляющей Да(?).
В научной литературе не удалось обнаружить объяснений причин возникновения систематической погрешности при обработке дискретного ЧМ сигнала конечной длительности методом ПК, в отличие от ПГ. Вероятно, это связано с тем, что метод ПК изначально является чисто эмпирическим и какой-либо строгой теоретической базы под собой, по сути, не имеет. Отсутствие аналитических выражений, дающих четкие численные соотношения между исходным сигналом и полученной огибающей кривой, - один из наиболее серьезных недостатков метода ПК, что признается и автором [6, 7].
Анализ результатов обработки ЛЧМ сигнала
В связи с необходимостью объяснения обнаруженных недостатков метода ПК представляется целесообразным провести анализ результатов обработки ЧМ сигналов, законы изменения фазы и частоты которых описываются аналитическими выражениями, например, сигналов с линейным законом частотной модуляции (ЛЧМ): ю (t) = ю j+ (ю 2 - ю j)t / T,
ф(?) = ф0 + ю it + (ю 2, t G [0; T], (11)
u(t) = sin ф^).
Результаты вычисления МЧ дискретных ЛЧМ сигналов, параметры которых приведены в табл. 1, с помощью ПГ и метода ПК представлены на рис. 4 и в табл. 2.
Для модельного ЛЧМ сигнала 1, как видно из табл. 2, погрешность вычисления МЧ методом ПК (0,3 % на отрезках [0; 0,2Т] и [0,8Т; Т], и 0,1 % на отрезке [0,2Т; 0,8Т]) оказалась значительно меньше аналогичной величины, полученной при использовании ПГ (9,5 и 0,9 % соответственно). Отметим, что, используя сглаживание значений МЧ, полученных с помощью ПГ, методом скользящего среднего или аппроксимируя данную зависимость линейной функцией по методу наименьших квадратов, можно уменьшить максимальную погрешность МЧ до 2 %, которая, однако, все еще на порядок превосходит ошибку метода ПК.
Для модельного ЛЧМ сигнала 2, как видно из табл. 2, погрешность вычисления МЧ методом
Таблица 1
Параметры модельных ЛчМ сигналов
Номер сигнала ш1, Гц ю2, Гц Т, с Число узлов временной сетки
1 10 20 1,0 32768
2 1000 3000 1,0 32768
Таблица 2
Относительная погрешность вычисления Мч сигналов
Метод вычисления МЧ Сигнал 1 Сигнал 2
Дю„ % Дю2, % Дю„ % Дю2, %
ПГ 9,3 0,9 251,7 1,1
ПК 0,3 0,1 11,7 8,5
Дш1 - максимальная погрешность значений МЧ по краям интервала анализа на отрезках [0; 0,2Т] и [0,8Т; Т], Дю2 - максимальная погрешность значений МЧ внутри интервала анализа на отрезке [0,2Т; 0,8Т].
ПГ на основном отрезке [0,2Т; 0,8Т], напротив, оказывается меньше (1,1 против 8,5 % метода ПК соответственно). При этом погрешность вычисления МЧ методом ПГ по краям интервала на отрезках [0; 0,2Т] и [0,8Т; Т] почти в два раза превосходит таковую у метода ПК (25,7 против 11,7 % соответственно). Отметим, что погрешность вычисления МЧ методом ПК увеличивается при увеличении частоты ЛЧМ сигнала.
Для объяснения обнаруженной особенности МЧ, вычисленной с помощью ПК, напомним, что в данном методе в качестве амплитудной составляющей сигнала а({) применяются результаты сплайн-интерполяции, для которой в качестве узлов интерполяции используют точки, соответствующие максимумам (минимумам) анализируемого сигнала. При этом априори предполагается,
что моменты времени, соответствующие достижению сигналом максимального (минимального) значений, известны абсолютно точно. Однако при использовании дискретной временной сетки ti = Дt • /, / = 0, N может оказаться, что текущее значение момента времени t. отличается от точного значения времени Т соответствующего данному максимуму (минимуму) анализируемого сигнала на величину 8t:
^ = Тт + (12)
тогда соответствующее значение дискретного сигнала, вычисляемого в соответствии с (11), и 1 Ф\, т. е. значения функции в узлах интерполяции оказываются заданными с некоторой погрешностью.
Получим оценку зависимости данной погрешности от частоты сигнала, разложив для это-
Рис. 4. Зависимость мгновенной частоты ЛЧМ сигнала от времени: а - модель 1; б - модель 2
10
10"
10"
10"
1tf
10
i 1 i i i / -j-i- \ ПК
- (ю,,-«Д Гц 2 I I 1 I 1 I Г '
1000
2000
4000
5000
300С
7000
Рис. 5. Зависимость погрешности вычисления МЧ от девиации частоты ЛЧМ сигнала
го функцию u (t) = sin ф(7) в окрестности точки Tm в ряд Тейлора:
u (Tm +St) = sinfaT +St)) =
= Sin^(Tm )) - sin(Ф(Tm ))
d 2ф ~dtT
St2
(13)
1 -®2
St2
Т 2
Из (13) видно, что при фиксированном значении Ы погрешность значений функции в узлах интерполяции (шум, обусловленный квантованием сигнала по времени) действительно увеличивается при увеличении девиации сигнала ю2 — ю1. Поэтому значения амплитудной составляющей сигнала а(г), вычисляемые с помощью сплайн-интерполяции в узлах временной сетки г. = Аг • /, / = 0, N, с неизбежностью не будут соответствовать априорным представлениям о том, что а(г1) = 1.
Приведенную выше аналитическую оценку также подтверждают результаты, представленные на рис. 5. (Здесь в качестве значения погрешности вычисления МЧ ЛЧМ сигнала принималась погрешность коэффициента прямой, аппроксимирующей табличную зависимость / = /(г.), значения которой находились для заданных значений девиации частоты ЛЧМ сигнала). Из рисунка видно, что при увеличении девиации частоты ЛЧМ сигнала погрешность МЧ, вычисленная в соответствии с ПГ, практически на меняется, в то время как аналогичная величина, вычисленная в соответствии с методом ПК, при увеличении девиации частоты (соответственно, шума сигнала, обусловленного квантованием по времени), напротив, увеличивается. Полученный результат позволяет предположить существование сильной
зависимости погрешности характеристик сигналов, вычисляемых методом ПК, от шумов, присутствующих в сигнале.
Для подтверждения данного предположения проведено статистическое моделирование, в котором модель ЛЧМ сигнала 2 дополнялась гауссовым шумом (г):
^ (г) = и (г) + ^ (г). (14)
При этом в процессе моделирования для каждого заданного значения отношения сигнал/шум: 1) генерировались 1000 независимых реализаций сигнала
»
(г1), I = 0, N, к = 1,1000, в соответствии с (14); 2) для каждой из реализаций и^ (¿.) были найдены табличные значения МЧ от времени — = (^); 3) вычислялись в соответствии с методом наименьших квадратов погрешности коэффициентов прямых Аак, аппроксимирующих табличные зависимости /^ = /^ (г.); 4) по массиву Аак вычислялось среднее значение погрешности углового коэффициента. Полученные зависимости погрешности МЧ от отношения сигнал/шум представлены на рис. 6.
Из рис. 6 видно, что, действительно, при низких отношениях сигнал/шум (< 22 дБ) погрешность нахождения МЧ методом ПК оказывается больше аналогичной величины, вычисляемой в соответствии с методом ПГ. При отношении сигнал/шум > 22 дБ, наоборот, погрешность нахождения МЧ методом ПК оказывается меньше аналогичной величины, вычисленной методом ПГ. Например, при отношении сигнал/шум > 40 дБ и более точность метода ПК выше, чем у метода ПГ, на постоянную величину ~5 %. Таким образом, полученные оценки погрешностей МЧ модельных сигналов следует принимать во внимание при выборе метода обработки реальных сигналов.
Рис. 6. Погрешность расчета линейного коэффициента ЛЧМ сигнала в зависимости от соотношения сигнал/шум
Проведенный сравнительный анализ результатов вычисления мгновенной частоты дискретных ЧМ сигналов конечной длительности методами преобразования Гильберта и прямой квадратуры позволяет сделать следующие выводы.
Погрешность вычисления мгновенной частоты модельных ЧМ сигналов методом прямой квадратуры оказывается на порядок меньше аналогичной величины, полученной методом преобразования Гильберта.
Каждому из изученных методов обработки сигналов присущи краевые эффекты, проявляющиеся в увеличении погрешности вычисляемых параметров сигналов на краях временного интервала.
Погрешность параметров ЛЧМ сигнала, вы-
числяемых в соответствии с методом прямой квадратуры, увеличивается при увеличении девиации частоты сигнала.
Метод прямой квадратуры оказывается более чувствительным к шумам, присутствующим в сигнале, чем метод преобразования Гильберта, т. е. намного сильнее зависит от соотношения сигнал/шум, чем ПГ.
Применение метода прямой квадратуры для анализа других модельных сигналов и нестационарных временных рядов, а также их сравнение с аналогичными результатами, получаемыми другими известными методами (например, вэйвлет-анализ, метод SSA и др.), являются предметом дальнейших публикаций.
СПИСОК ЛИТЕРАТУРЫ
1. Gabor, D. Theory of communication [Текст] / D. Gabor // J. of the Institution of Electrical Engineers. -1946. -Vol. 93 (III) -P. 429-457.
2. Рытов, С.М. Введение в статистическую радиофизику [Текст] / С.М. Рытов. -М.: Наука, 1966. -С. 452.
3. Вакман, Д.Е. Разделение частот в теории колебаний и волн [Текст] / Д.Е. Вакман, Л.А. Вайнштейн. -М.: Наука, 1983. -С. 288.
4. Поршнев, С.В. Радиолокационные методы измерений экспериментальной баллистики [Текст] / С.В. Поршнев.- Екатеринбург: УрО РАН, 1996.- С. 212.
5. Поршнев, С.В. Физическое содержание понятий «огибающая» и «мгновенная частота широкополосного аналитического сигнала» [Текст] / С.В. Поршнев // Электромагнитные волны и электронные системы. -2001. -Т. 6. -№ 1. -С. 48-55.
6. Huang, N.E. The Hilbert-Huang transform and its applications [Текст] / N.E. Huang; ed. Samuel S.P. Shen. -World Scientific Publishing Company Pte. Ltd. 5 Tuck. Link, Singapore, 2005. -P. 324.
7. Huang, N.E. On instantaneous frequency [Текст] / N.E. Huang [et al.] // Advances in Adaptive Data Analysis. -World Scientific Publishing Company, 2009. -Vol. 1. -№ 2. -P. 177-229.
8. Коржик, В.И. Огибающая сигнала и некоторые ее свойства [Текст] / В.И. Коржик // Радиотехника. -1968. -Т. 23. -№ 4. -С. 1-6.
9. Коржик, В.И. Расширенное преобразование Гильберта и его применение в теории сигналов [Текст] / В.И. Коржик // Проблемы передачи информации. -1969. -Т. 5. -Вып. 4. -С. 3.
10. Рандал, Р.Б. Частотный анализ [Текст] / Р.Б. Рандал. -ДК-2600 Глоструп, Дания: А/О К. Ларсен и сын, 1989. -С. 389.
11. Shekel, G. On the Term Instantaneous Frequency [Текст] / G. Shekel // Proc. I.R.E. -1954. -Vol. 42. -№ 6. -P. 1024.
12. Тихонов, В.И. Один способ определения огибающей квазигармонических флюктуаций [Текст] / В.И. Тихонов // Радиотехника и электроника. -1957. -Т. 2. -№ 4. -С. 562.
