Модификация преобразования Вигнера-Виля для анализа интерферометрических данных газодинамических процессов Текст научной статьи по специальности «Компьютерные и информационные науки»

Радиофизические измерения Вестник Нижегородского университета им. Н.И. Лобачевского, 2011, № 5 (3), с. 95-103

УДК 621.396

МОДИФИКАЦИЯ ПРЕОБРАЗОВАНИЯ ВИГНЕРА-ВИЛЯ ДЛЯ АНАЛИЗА ИНТЕРФЕРОМЕТРИЧЕСКИХ ДАННЫХ ГАЗОДИНАМИЧЕСКИХ ПРОЦЕССОВ

© 2011 г. С.Ю. Лупов, В.И. Кривошеев

Нижегородский госуниверситет им. Н.И. Лобачевского

1ироу@г£ ипп. ги

Поступила в редакцию 12.05.2011

Приведен алгоритм, позволяющий оценить частотно-временное распределение энергии (ЧВРЭ) коротких широкополосных частотно-модулированных сигналов. Приведен пример вычисления ЧВРЭ интерферограммы, полученной в эксперименте по метанию металлической пластины продуктами взрыва заряда из тротила.

Ключевые слова: преобразование Вигнера-Виля, частотно-временное распределение, газодинамический процесс, интерферограмма.

Введение

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

Отметим, что задача непрерывного измерения параметров движения снаряда в начальный момент выстрела до настоящего времени не решена [1], т.к. кинематические характеристики и, прежде всего, скорость на начальном этапе выстрела меняется в широком диапазоне значений за короткий временной промежуток.

Для решения таких задач в конце XX века получил широкое распространение радиоин-терферометрический метод измерений. Предложенные в [1, 2] алгоритмы для анализа ин-терферометрических данных не всегда эффективны, т.к. не учитывают многомодовость сигнала, возникающую из-за переотражений электромагнитной волны между движущимся объектом и антенной радиоинтерферометра. Для анализа таких сигналов могут быть использованы различные методы частотно-временного анализа [3, 4], в том числе и вэйвлет-анализ [5], но, т.к. в основе таких методов, в явном или неявном виде, лежит дискретное преобразование Фурье, точность анализа зависит от количества периодов на исследуемом участке интер-ферограммы. В эксперименте по метанию

стальной пластины продуктами взрыва заряда взрывчатого вещества [6] интерферограмма имеет около 40 периодов, но наибольший интерес исследователей представляет начальный участок интерферограммы, на котором мгновенная частота основной компоненты меняется от нулевой до установившейся примерно за 10 периодов. Изменение частоты происходит ступеньками [2, 6]. В данной статье приводится алгоритм, позволяющий оценить ЧВРЭ сигнала на данном участке.

1. Описание алгоритма

Псевдопреобразование Вигнера

В основу алгоритма положено преобразование Вигнера-Виля, описанное в [3]:

/ ) = 15(х + У2 )'(х“ У2 )ехр(- ]2тф )*.

—ад

Функция ) называется ЧВРЭ сигнала, может принимать только действительные значения (включая отрицательные).

Несмотря на высокое разрешение, как по частоте, так и по времени, в распределении могут образовываться побочные частотные компоненты, мешающие анализу сигнала [3, 7]. Это связано с нелинейностью преобразования.

Существует несколько методов, позволяющих уменьшить интенсивность побочных компонент, используя определенные процедуры усреднения. Один из них - использование окна Ь(р) во временной области. В результате получается так называемое псевдопреобразование Вигнера (ППВ) [4]:

то при t0 ^ ад ППВ переходит в обычное преобразование Вигнера-Виля. При уменьшении t0 интенсивность побочных продуктов ослабляется, платой за это является ухудшение частотного разрешения.

При анализе оцифрованного сигнала ППВ удобнее вычислять с помощью быстрого преобразования Фурье (БПФ) в скользящем окне. Для этого перед выполнением процедуры БПФ отрезок сигнала s[n], выделенный скользящим окном с размером M отсчетов, преобразуют по следующему алгоритму.

1) если размер окна нечетный, то

S[n]= s[n\s*[M — n — 1], n = 0,1,...,M — 1,

для четного размера окна

5[п] = 5[и]' 5*\М -п - 2], п = 0,1,...М - 2, 5 М -1] = s\M -1] • 5 * [м -1}

2) чтобы результат процедуры БПФ получился действительным, необходимо перед вычислением БПФ выполнить циклическую перестановку полученного сигнала Sl[n] влево на (М-1)/2 отсчетов (если М - нечетное) или на М/2-1 (если М - четное).

При построении ЧВРЭ все значения на шкале частот следует разделить на 2.

В качестве примера построим ЧВРЭ для смоделированного сигнала:

^[п] = ехр(/(2тс • 0.05п + 1008Іп(2тс • 0.0005?))) + + ехр(/ (2 л • 0.1п + 2008Іп(2тс • 0.0005?))),

N = 0, 1, ..., 3999,

который состоит из двух ЧМ компонент, мгновенная частота одной из них меняется по синусоидальному закону в диапазоне от 0 до 0.1, а другой - от 0 до 0.2.

Рис. 1. ЧВРЭ, вычисленное с помощью псевдопреобразования Вигнера (M=500)

Рис. 2. Фурье-спектрограмма (M=500)

Pra. 3. ЧBPЭ, вычисленное с помощью псевдопреобразования Вигнера (M=2000)

Pra. 4. ЧBPЭ, вычисленное с помощью псевдопреобразования Вигнера (M=500)

На рис. 1 представлено ЧВРЭ, полученное с помощью ППВ с размером окна M=500 отсчетов. По оси абсцисс отложено время, по оси ординат - частота. Более темные участки распределения соответствуют большей интенсивности.

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

Качественно можно видеть, что ЧВРЭ ППВ (рис. 1) имеет более высокое частотно-временное разрешение по сравнению со спектрограммой (рис. 2). При увеличении размера окна количество и интенсивность побочных компонент в ЧВРЭ ППВ увеличиваются, что может затруднить выделение основных компонент (рис. 3).

Рассмотрим другой пример.

Смоделированный сигнал:

s\n\ = ехру(2тс • 0.05n + sin(2rc • 0.0]n)))+

+ exp(j(2^ • 0.1n + 2 • sin(2ft- 0.0]л))),

N = 0, 1, ..., 999 (1)

состоит из двух ЧМ компонент, мгновенная частота одной из них меняется по синусоидальному закону в диапазоне от 0.04 до 0.06, а другой - от 0.08 до 0.12, в отличии от предыдущего сигнала диапазон частот каждой ЧМ компоненты уменьшился в 5 раз, а период закона изменения частоты - в 20 раз.

На рис. 4, 5 представлены ЧВРЭ, полученные с помощью ППВ с различными размером окон. При размере окна в 500 отсчётов (рис. 4) количество и интенсивность побочных компонент в ЧВРЭ затрудняют выделение ЧМ компонент, присутствующих в сигнале. Уменьшение размера окна до 50 отсчетов (рис. 5) позволяет оценить структуру сигнала, но разрешение по частоте может быть недостаточно для анализа сигнала.

Модификация псевдопреобразования Вигнера

Чтобы исследовать сигналы, подобные (1), была сделана модификация алгоритма вычис-

Рис. 5. ЧВРЭ, вычисленное с помощью псевдопреобразования Вигнера (М=50)

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

1) Отрезок сигнала s\n], выделенный скользящим окном размером М отсчетов, преобразуем по следующему алгоритму:

^[п] = ь\М1 + п]-я*[М1 -п], п = 0,1,...,М1, (2)

М-1 М

где М =-------- - для нечетного М, М =--------1

1 2 1 2

- для четного М.

2) Используя алгоритм оптимизации, аппроксимируем полученный сигнал .^[п] линейной комбинацией р комплексных экспонент (не обязательно ортогональных на данном отрезке) с различными частотами и амплитудами Ак, так, чтобы относительная ошибка аппроксимации между моделью и сигналом была минимальной.

Относительная ошибка аппроксимации вычисляется по формуле:

1

XI 5і["]_

8 =

п=0

1

XIФ]

(3)

г

где ~1[п] = ХАк

п=0

к=1

и 1_ 1

X АкС1к = X (^е(^1 [и]) со82/и + іш(.ї[и])8Іп2/1и)

к=1 п=0

, (4)

Р М1 , ч

X Аксрк = ХМ51[и])со827/ри + іш(і[и])8Іп2л/ріи)

Ак

_к=1 п=0

где

С

тк

‘■-Ч

= Х “Ф^т - /к )«) =

/т * /к;

Г 1 , 8ш[(М 1 + 0.5)2к(/т - /к )]

модель полученного

сигнала, / - параметр поиска, который может принимать действительные значения от -0.5 до 0.5.

Оптимальное значение Ак находится из системы линейных уравнений:

2 ^тК/т - /к ))

М + 1, /т = /к .

3) Полученные значения частот делим на 2.

4) Значения амплитуд отображаются на плоскости время-частота. По оси абсцисс откладывается номер отсчета сигнала, соответствующий центру скользящего окна, по оси ординат - линейная цифровая частота, по оси аппликат на частотах /к отображаются значения амплитуд Ак оттенками серого.

Блок-схема алгоритма

На рис. 6 представлена блок-схема алгоритма.

Отрезок сигнала s\n], выделенный скользящим окном с размером М отсчетов (блок 1), преобразуют в 51 [п] (блок 2).

Дальнейшая обработка преобразованного сигнала происходит в цикле (блок 3). С помощью генератора случайных чисел (блок 4) задаются начальные значения всех цифровых частот/к в диапазоне от -0.5 до 0.5. В блоке 5 проводится поиск оптимальных параметров при которых относительная ошибка аппроксимации, вычисляемая в целевой функции (блок 6), имеет

2

2

,г5[0]^.5[М-1] Преобразование сигнала ^[и]^51[п], используя соотношение (2)

2

•я[0],..., 51 [М-1]

м]

Начальные случайные значения /ь /2,., /р 4 ЦЕЛЕВАЯ ФУНКЦИЯ 1. Вычисляет значения Ак, решая систему линейных уравнений (4)

і г /1,.,,/р 2. Вычисляет относительную

Поиск оптимальных ошибку аппроксимации е

параметров (симплексный алгоритм) 5 8 по формуле (3) 6

51 [0],... 5і[М-11

ш Цикл выполняется заданное число раз 3 V

минимальное значение. В процессе поиска значения цифровых частот могут перейти за границы заданного диапазона, поэтому полученные данные приводятся в диапазон от -0.5 до 0.5. Весь этот процесс происходит многократно, сохраняя значения поисковых параметров, при которых целевая функция имеет наименьшее значение.

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

Особенности реализации алгоритма

Как правило, в сигналах присутствует постоянная составляющая. Чтобы ее учесть, в экспо-

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

Для поиска оптимальных параметров хорошо зарекомендовала подпрограмма Е04ССЕ [9] библиотеки КЛвЫБ, в которой реализован симплексный [10] алгоритм оптимизации.

Особенности вычисления частоты и амплитуды

При вычислении ЧВРЭ широкополосных ЧМ сигналов алгоритмами, основанными на оценке параметров синусоидальной или экспоненциальной модели методом оптимизации,

0.15-

0.10- \ А. А А А І-. .•. Г\. • ■< •' < • < • -І • ■ • ‘ . V • ‘ ‘ V '■ " ' . » • •. ‘ \ . и ‘ .* •. У’ \ і \ ■ ■ А А ■! '• 1 " і •V . < V . V . * . V . ' / . * /•* =• і ’ ' Л'

! і /\\ ■■■■ '-А ■ • '■ \у і • / '• : • '•'. * •, '•'. ■ • ;', ‘ - . ' и /'А •*. / .*• ■ ?' : »• ; '• .*

0 05- '■■■ 4 .* •**'.' ' /***•• 4 / ’ Л • >. ■ - >■ ■ і. > ■ /!• Л •\ \ \ \ Ч/' Ч-: V* •\ •• ■ ■ : - Л'. К;'. ••/’ • V- ■;■/■• V• чД’ д* > , **».• < / > • /■’ } і • Ь Л • • :.Г •V *. -V •. V

0.00-

100 200 300 400 500 Затріев 600 700 800 900

Рис. 7. ЧВРЭ, вычисленное с помощью модифицированного преобразования Вигнера-Виля (М=15, Р=4)

0.15-

0.101 < Л Л :Ч . ..*-4 л’ *' ; л' > л ► л V л* > л’ N Л' ' Л* N л* “ 'Л’' ^ * Лч^ Л'о. Лч'- Гл« Г*с Лч'- Л. Л. Л: =■ ► л’ а Л > Л А V ■>' ' ■> N л* • •* • • І •• і О. ІЇС- Лїо

її: 0.05* .М. ./*Л. А. .А. /Гл. .а. .А ,/л. а. /?' -я /V ^ ^ /’* ^ /• ^ X Ъ Р Ь р \ ;? ^ л? і ^ ^ ^ ^ ^ •<* ^ ^ ^ &

0.00-

100 200 300 400 500 600 Батріев 700 800 900

Рис. 8. ЧВРЭ, вычисленное методом оценки оптимальных параметров экспоненциальной модели в скользящем окне (М=15, Р=3)

возникает побочный эффект, когда амплитуда двух или более частотных компонент во много раз превышает амплитуду сигнала. Если ЧМ компонента определена на коротком отрезке, то ее можно аппроксимировать с высокой точностью несколькими комплексными экспонентами с близкими частотами и амплитудами, во много раз превышающими амплитуду исходного сигнала, сумма начальных фаз этих компонент близка к п. При сложении этих компонент мгновенное значение частоты и амплитуды на заданном отрезке меняется по закону [11]:

1 х[і ]• ]- у[і] х[і ]

2п

’[; ]

+у2 И

(5)

Л2[] = х2[]+у2[], (6)

где

х[]=Е Лк со Фп/ і+Фк),

к

уИ=Е Лк *іп(2п/кі + Фк ),

к

х'[]= -2пЕЛк/к ^п(2п/к і +Фк),

к

у ’[і]=2пЕ Лк/ксо і2/і+Фк) •

к

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

х

□ 500 1000 1500 2000 2500 3000

ватрЬб

Рис. 9. Квадратурный сигнал, полученный на выходе интерферометра

0.020-

0.005- ;

о.ооо-|- _ _ ( _ ........................., , , ,

0 500 1000 ' ' ' 1500 ’ ' 2000 ' ' 2500

8атр!е8

Рис. 10. ЧВРЭ, вычисленное с помощью модифицированного преобразования Вигнера-Виля интерферограммы (М=200, Р=2)

шает заданное значение1, сумма фаз близка к п, то следует объединить эти компоненты в одну, пересчитав амплитуды и частоты по формулам (5) и (6). Для алгоритма, приведенного в [8], значение t равно половине размера окна (не округляя до целого). Для алгоритма, приведенного в данной статье, значение t = 0. В этом случае соотношения (5) и (6) можно переписать:

, X ^кАк

1 [°]=2П ХТ'

к

А[°] = Х Ак .

к

1 В реализованной программе для алгоритмов приведенных в [8] и в данной статье это значение равно

0.0001 (подбиралось экспериментально).

Выбор количества частотных компонент в модели

Количество частотных компонент в модели выбирается, исходя из априорной информации о сигнале.

2. Тестирование алгоритма

В качестве примера на рис. 7 представлено ЧВРЭ смоделированного на компьютере сигнала (1), вычисленное с помощью модифицированного преобразования Вигнера-Виля (МПВВ) (размер окна М=15 отсчетов, количество частотных компонент, задаваемых для поиска, Р=4). Несмотря на малый размер окна, удалось не только разделить частотные компоненты, но и достаточно точно описать закон изменения частоты каждой из них.

Для сравнения на рис. 8 представлено ЧВРЭ, вычисленное с помощью метода оценки опти-

мальных параметров экспоненциальной модели (ООПЭМ) в скользящем окне, приведенного в [8], на котором хорошо различимы частотные компоненты, но, в отличии от МПВВ, на них заметны ложные ступеньки.

3. Применение модифицированного алгоритма Вигнера-Виля для анализа экспериментальных данных

Рассмотрим пример анализа данных, полученных в опыте по регистрации процесса метания стальной пластины продуктами взрыва заряда из тротила [2, 6].

На рис. 9 приведена запись выходных сигналов интерферометра, на рис. 10 - ЧВРЭ, вычисленное с помощью МПВВ (размер окна М = 200 отсчетов, количество частотных компонент в модели Р = 2). На ЧВРЭ отчетливо видна частотная компонента, соответствующая скорости пластины. Соотношение, приведенное в [2]

V = -/

2

позволяет пересчитать мгновенную частоту в скорость V пластины. В данном соотношении / - реальная доплеровская частота, которая вычисляется умножением цифровой частоты на частоту дискретизации, - - длина волны излучаемой радиоинтерферометром.

В ЧВРЭ отчетливо видны ступеньки на начальной стадии движения пластины (первые 1000 отсчетов), которые действительно присутствуют в сигнале [2].

Выводы

Описанный алгоритм предназначен для вычисления ЧВРЭ широкополосного ЧМ сигнала.

Алгоритм основан на вычислении псевдопреобразования Вигнера, в котором заменили БПФ на оценку оптимальных параметров экспоненциальной модели методом оптимизации, что позволило анализировать ЧМ сигналы с небольшим количеством периодов. Малый размер скользящего окна увеличил разрешение по времени, при этом разрешение по частоте осталось высоким.

Анализ интерферограммы, полученной в опыте по регистрации процесса метания стальной пластины продуктами взрыва заряда взрывчатого вещества позволил оценить скорость пластины на начальной стадии движения.

Недостатком алгоритма является длительное время вычисления.

Список литературы

1. Поршнев С.В. Радиолокационные методы измерений экспериментальной баллистики. Екатеринбург: УрО РАН, 1999. 211 с.

2. Канаков В.А., Лупов С.Ю., Орехов Ю.И. Родионов А.В. Методы извлечения информации о перемещении границ раздела в газодинамических экспериментах с использованием радиоинтерферометров миллиметрового диапазона длин волн // Изв. вузов. Радиофизика. 2008. Т. 51, № 3. C. 234-246.

3. Коэн Л. Время-частотные распределения: Обзор // ТИИЭР. 1989. Т. 77, № 10. С. 72.

4. Лазоренко О.В., Черногор Л.Ф. Сверхширо-кополосные сигналы и физические процессы. 2. Методы анализа и применение // Радиофизика и радиоастрономия. 2008. Т. 13, № 4. C. 270-322.

5. Поршнев С.В. Применение непрерывного вейвлет-преобразования для обработки широкополосных частотномодулированных сигналов // Вычислительные методы и программирование. 2003. Т. 4, № 3. С. 104-116.

6. Родионов А.В., Канаков В.А., Лупов С.Ю. Методы обработки результатов радиоинтерферомет-рических измерений параметров газодинамических процессов // Экстремальные состояния вещества. Детонация. Ударные волны: Труды Междунар. конф. «7 Харитоновские тематические научные чтения». Саров, 2005. С. 680-685.

7. Шкелев Е.И., Кисляков А.Г., Лупов С.Ю. Методы ослабления эффектов интермодуляции в распределении Вигнера-Вилля // Изв. вузов. Радиофизика. 2002. Т. 45, № 5. C. 433-442.

8. Лупов С.Ю., Серебряков А.М., Фрадкина Е.П. Оценка оптимальных параметров экспоненциальной и синусоидальной моделей отрезка дискретного сигнала // Вестник ННГУ. 2011. № 2. С. 60-68.

9. http://www.nag.co.uk/numeric/Fl/manuai/pdf/E0 4/e04ccf.pdf (дата обращения: 15.05.2011)

10. Nelder J.A., Mead R. A simplex method for function minimization // Computer Journal. 1965. V. 7. P. 308-313.

11. Гоноровский И.С. Радиотехнические цепи и сигналы: Учебник для вузов. М: Радио и связь, 1986. 512 с.

A MODIFIED WIGNER-VILLE TRANSFORM TO ANALYZE INTERFEROMETRIC DATA

OF GAS-DYNAMIC PROCESSES

S.Yu. Lupov, V.I. Krivosheev

An algorithm to estimate time-frequency energy distribution (TFED) of short broadband frequency-modulated signals is presented. We give an example of TFED interferogram obtained in the experiment on throwing a metal plate by TNT explosion products.

Keywords: Wigner-Ville transform, time-frequency distribution, gas-dynamic process, interferogram.