CC BY
57
Радиофизика
Вестник Нижегородского университета им. Н.И. Лобачевского, 2011, № 2 (1), с. 71-80
УДК 621.391.037.372
ОЦЕНКА ОПТИМАЛЬНЫХ ПАРАМЕТРОВ ЭКСПОНЕНЦИАЛЬНОЙ И СИНУСОИДАЛЬНОЙ МОДЕЛЕЙ ОТРЕЗКА ДИСКРЕТНОГО СИГНАЛА
© 2011 г. С.Ю. Лупов, А.М. Серебряков, Е.П. Фрадкина
Нижегородский госуниверситет им. Н.И. Лобачевского
lupov@rf. unn. т
Поступила в редакцию 20.05.2010
Приведен алгоритм оценки оптимальных параметров синусоидальной и экспоненциальной моделей дискретного сигнала. Показана возможность использования данного алгоритма для оценки частотновременного распределения квазигармонического сигнала. Приведен пример вычисления частотновременного распределения реального сигнала, полученного на установке «Маятник».
Ключевые слова: оценка спектра, частотно-временное распределение, оптимизация.
Введение
Одной из актуальных задач цифровой обработки данных является задача определения параметров зарегистрированных сигналов [1]. Наибольший интерес представляют не единичные значения параметров, а их функциональные зависимости от времени. К этому классу задач относится вычисление частотно-временного распределения [2, 3].
Несмотря на множество различных алгоритмов, до сих пор актуальной остается задача вычисления частотно-временного распределения сигналов со слабой частотной модуляцией на коротком временном интервале [4].
В данной статье:
• приведен алгоритм оценки оптимальных параметров синусоидальной и экспоненциальной моделей дискретного сигнала методом оптимизации;
• описана методика вычисления данным алгоритмом в скользящем окне частотно-временного распределения квазигармонического сигнала;
• приведены результаты тестирования алгоритма на сигналах, полученных с помощью компьютерного моделирования и на лабораторной установке «Маятник» [5].
1. Описание алгоритма
За основу был взят алгоритм, описанный в [6, 7].
Для оценки спектра будем аппроксимировать отрезок дискретного сигнала $ линейной комбинацией Р синусоид или комплексных экспонент (не обязательно ортогональных на дан-
ном отрезке) с различными частотами гак, амплитудами Ак и фазами фк так, чтобы относительная ошибка аппроксимации между моделью и сигналом была минимальной.
Относительная ошибка аппроксимации вычисляется по формуле:
в =
N
X ^ ~
і=1_________
N , |2
і=1
(1)
где
Si = X А С08(шіі + ф£ ) - синусоидальная
к=1
модель действительного сигнала,
(2)
р
Л^кі+ф к)
- экспоненциальная
~ = X Ае'
к=1
модель комплексного сигнала. (3)
Оптимальные значения Ак и фк можно назвать амплитудами и фазами спектра на частотах лишь условно, т.к. данная модель учитывает особенности сигнала только с дискретным спектром, когда Р больше числа частотных компонент в сигнале или равно ему. Однако, как будет показано ниже, использование данного алгоритма может помочь оценить динамику отдельных частотных компонент в квазигармоническом ЧМ-сигнале.
Минимизировать относительную ошибку аппроксимации (1), зависящую от 3Р параметров Ш£, Ак и фк, аналитически невозможно, поэтому для поиска оптимальных параметров модели будем использовать алгоритм оптимизации.
2
Легко показать, что значения амплитуд Аь и фаз фь можно вычислить аналитически, что сокращает количество поисковых параметров в 3 раза.
Вычисление значений амплитуд и фаз в модели действительного сигнала. Для действительного сигнала значения частот положительны.
Подставим выражение (2) в (1):
N-1/
8 =
i г О
і S, -і Ak cos(qIі + фк )
A2
k=1
N-1
і S2
i=0
и, проделав несложные тригонометрические преобразования, получим
N-lf P \2
I Si -!(Xi • cos^ii - Yk sin&ki)
і=о v i=l
8 =
P Ex, k=1 cc1i P "lYt k=1 N -1 cslk = IS i=0 cos Q1i
Xk Jl ccPk P -'Ey, k=1 N -1 • csPk = IS, i=0 • cos Qpi
NTM P kX .. k sc1k P -EYi k=1 N -1 sslk = IS • i=0 sin Q1i
P 'LXi k=1 ScPk P -lYi k=1 1 Si = .. k ... P Co Co • Sin QPi
где
cc
mk
N-1 = і:(<
cOSQ^i • cOS^ii
m
0=
i=0
1 , sin((N - 0 . 5)-(Qm -Qi))
+
+
2 4 •sin((Qm -Qi V2)
sin ((N - 0.5)-(ram + юк))
+-----------„ ’ v m, . e”, Qm *Qk
4 •sin((Qm + Qi V2)
N +1 sin Qm - sin((2N -1) • Qm)
4 • sin
Qm = Qi * 0
Qm = Qi = 0
N-1
csmi = I(cOSQm • Sin Qii)г i=0
c°s((cam -Qi)/2)-cos((N -0.5>(ют -Qi)) +
4 •sin((Qm -Qk V2)
+ cos((Qm + Qk V2)- cos((N - 0. 5MQ m + Qk)) ’ 4 •sin((Qm +QiV2)
Qm *Qk
cosrnm - cos((2N - l)^)
4 • Sin (Qm )
0,
, Qm =Qi * 0 Qm = Qk = 0
(4)
sc
N-1 mk = I(si
Sin Qmi • cOSQki) = cS
km
N-1
IS2
2=0
где Xk = Ak cos9k и Yk = Ak sin ф* .
Для нахождения оптимальных параметров Xk
ds ds
и Yk вычислим производные ---------, ----- (m =
dXm dYm
= 1, 2, ... , P) и приравняем их к нулю, получится система 2P линейных уравнений:
i=0 N-1
SSmk = I(sinQmi •SinQki)г
i=0
sin((N - 0.5)^(rom -Qk))
4 •sin((Qm -Qk V2)
sin((N - °.5MQm + Qk)) ................
а ~Л , Qm *Qk
4 • sin((Qm + QkУ2)
N +1 sin Qm - sin((2N - !)• Qm)
2 4 • sin(rom)
Qm = Qk * 0
, (5)
0, ют =ак = 0
из которой найдем Хь и Уь, такие, чтобы относительная ошибка аппроксимации е была минимальной для заданных частот юь. Оптимальные значения амплитуд Аь и фаз фь будут получены из выражений
Фк = Агсіе-1*
Ai =,1 Xl + Yi2
X
Вычисление значений амплитуд и фаз в модели комплексного сигнала. Подставим выраже ние (3) в (1), получим
р / Л( р
8 г -
N-1
і
= 0
к=1
s, - XAkej(fflki+%) її S* - pAke~J(икі+фк)
к =1
N-1
і |Si|2 i=0
Сделав замену hi = Aie]i®k , получим
8 = -
N-1
I
i=0
P
I
к=1
к = Ake
\ f
Sі -і hj
P
<* X—' » *
S* - іh*e Q
к=1
N -1
(б)
і |Si|-
і г 0
P
2
Для нахождения оптимальных параметров Нк вычислим производные ------- (т = 1, 2, ..., Р) и
приравняем их к нулю, получится система Р линейных уравнений:
N-1
hlgl1 + h2gl2 + ... + hpglP = I Si
~JQll
i=0 N-1
h1g21 + h2g22 + ... + hpg2P = I Sie У 2
i=0
N-1
hlgP1 + h2gP2 + ... + hpgpp = I Sie J“Pl
i=0
где
N -1
gmk I e
j(qi-Qm ) _
i=0
1 - e
j(“k-Qm )N
(“i
Qk *Qm
N -1
cmk = I cOs((“m +Qk )')г
i=0
1 , sin((N - 0.5)(“m +Qk))
— + 2 N
2 2sin ((“m +Qk V2)
Qm +Qk * 0 ,
Qm +Qk = 0
mk
IN -1
= I Sin ((“m +Qk )•1 )г
,(7)
cO s((“ m +Qk V2)-cOs((N - 0.5 XQ m +Qk ))
2 • Sin ((Qm +Qk V2)
Qm +Qk * 0
Qm +Qk = 0.
Оптимальные значения амплитуд Аь и фаз фк будут получены из выражений
}к ют
из которой найдем значения Ьк, такие, чтобы ошибка аппроксимации е была минимальной для заданных частот юк.
Для реализации алгоритма на компьютере выражение (6) удобнее представить в тригонометрическом виде. Формула (6) примет вид:
. Imhk Yk
фі = Arctg—-к- = Arctg-f-.
k Re h X
Блок-схема алгоритма. На рис. 1 представлена блок-схема работы алгоритма.
С помощью генератора случайных чисел (блок 1) задаются начальные значения всех
N -1
I
i =0
8 = ■
(
V
(
Re si -I(Xk cos“ki - Yk sin “ki) + ImSt - I(Xi sin q,і + Y, cosq,і)
k = 1 / v k = 1
V
N -1
(8)
I Si Ґ
i г 0
где Xk = Re hk = Ak cos фk и Xk = Im hk = цифровых частот ak в диапазоне от 0 до ж для
= Ak sin ф^
мет вид:
P P N-1
IXlc1l - IYlS1l = I(ReSi cOSQ1i - ImSi SinQ1l) l=1 l=1 i=0
действительного или от -тс до тс для комплекс-Тогда система линейных уравнений (7) при- ного сигнала. В блоке 2 производится поиск
оптимальных параметров ак, при которых относительная ошибка аппроксимации, вычисляемая в целевой функции (блок 3), имеет минимальное значение. В блоке 4 реализована визуализация исследуемого и аппроксимированного сигналов (в виде графиков), найденных значений частот, амплитуд и фаз (в виде таблицы) и т.д. В процессе поиска значения цифровых частот могут перейти за границы заданного диапазона, поэтому полученные данные при отображении приводятся в диапазон от 0 до тс для действительного или от -тс до тс для комплексного сигнала.
Весь этот процесс происходит многократно, сохраняя значения поисковых параметров, при которых целевая функция имеет наименьшее где значение.
P P N -1
IXlcPl - IYkSPk = I(Re Si cosQp - ImSi sin QPi) l=1 l=1 i=0
, (9)
N-1
IXkslk + IYicli = I(ReS1 smQll +ImSi cosQll) l=1 l=1 i=0
P P N -1
IXkSPk + IYicPi = I(ReSi SinQPl +ImSi coSQPi') l=1 l=1 i=0
— <
e
i=0
— <
0
m
p
p
Цикл остановится, если выполнился заданное число раз, ц| е стала меньше заданного значения или нажата кнопка «Stop» q
Рис. 1. Блок-схема алгоритма
Особенности реализации целевой функции. При вычислении целевой функции предварительно сравниваются друг с другом входные параметры. При совпадении двух или более параметров системы линейных уравнений (5) и (7) становятся вырожденными. Чтобы избежать этого, находим параметры с равными значениями (с точностью 10-12), оставляя только один и сокращая количество уравнений в системе.
Дальнейшее изменение целевой функции связано с особенностями исследуемых сигналов. Как правило, в сигналах присутствует постоянная составляющая. Чтобы ее учесть, в целевой функции добавляется параметр raP+i=0, который не является поисковым. В комплексных сигналах реальная и мнимая части редко бывают строго в квадратуре, поэтому добавляется еще P параметров со значениями fflk+P+i = = -fflk-
Введение дополнительных частотных компонент в целевой функции позволяет сократить количество поисковых параметров, что значительно сокращает время поиска.
Для поиска оптимальных параметров в блоке 2 (рис. 1) использовались:
• подпрограмма E04CCF [8] библиотеки NAGLIB;
• виртуальный прибор LabVIEW Downhill Simplex nD, модифицированный для работы со сложными функциями [9];
• подпрограмма E04JAF [10] библиотеки
клеив.
В первых двух реализован симплексный [11], в третьей - квазиньютоновский [12] методы оптимизации. Наименьшее время поиска оптимальных значений параметров, при которых относительная ошибка аппроксимации (8) имеет глобальный минимум, показала программа, использующая E04CCF. Чуть медленнее работает программа, написанная на языке программирования ЬаЬУІЕ"^ При использовании подпрограммы E04JAF время поиска глобального минимума целевой функции в среднем увеличивается в несколько раз.
Выбор количества частотных компонент в модели. Количество частотных компонент в модели выбирается исходя из априорной информации о сигнале.
3. Тестирование алгоритма
Рассмотрим пример анализа тестового комплексного сигнала, записанного выражением:
( 4 Л
X cos(2/ )+p1
V k=1
+
Г 4 ^
+ j Xsin(2/) + P2
V k=1
i = 0, 1,...,399
Рис. 2. Результат аппроксимации сигнала без шума экспоненциальной моделью
где Д = 0.001, f2 = 0.01, Дз = 0.011, Д = 0.022, р: и р2 - два независимых источника белого гауссова шума. Чтобы проверить разрешающую способность алгоритма по частоте, значения частот второй и третьей компонент выбраны так, чтобы период биений, создаваемый ими, в 2.5 раза превышал длительность сигнала. Для проверки способности алгоритма оценивать сверхнизкие частоты период первой компоненты в 2.5 раза превышает длительность сигнала.
На рис. 2 представлены результаты аппроксимации сигнала без шума. Количество поисковых параметров равно 10 (в целевой функции задавались еще 11 непоисковых параметров с нулевой и противоположными по знаку частотами). Относительная ошибка аппроксимации (1) получилась равной 6.8-10-16. Время поиска не превышает 1 минуты1. На верхнем графике изображена реальная, на нижнем - мнимая части исходно-
1 Тестирование проводилось на ноутбуке 8аш8ищ МР-К560-Б802Яи при использовании подпрограммы поиска оптимальных параметров Е04ССР.
го и аппроксимированного сигналов (аппроксимированный и исходный сигналы полностью совпадают). В таблице справа представлены вычисленные параметры (в колонке А - амплитуды, в колонке Е - цифровые линейные частоты, в колонке Е - фазы). Значения 2-й, 4-й, 5-й и 6-й строк полностью совпадают с параметрами исходного сигнала.
На рис. 3 представлены результаты аппроксимации сигнала с дисперсией шума, равной 1. Количество поисковых параметров
- 10. Относительная ошибка аппроксимации (1) получилась равной 0.29. Среднее время поиска увеличилось до 30 минут. На графиках исходный сигнал изображен более светлым оттенком, аппроксимированный - более тёмным. Значения 2-й - 5-й строк в таблице (рис. 3) соответствуют 1 -й - 4-й частотным компонентам исходного сигнала. Ошибки по частоте 2-й, 3-й и 4-й компонент сигнала не превышают 0.0001 (или 1%). Ошибка по частоте низкочастотной компоненты не превышает 0.0002 (или 20%). Ошибка по амплитуде не превысила 13%.
Рис. 3. Результат аппроксимации сигнала с шумом экспоненциальной моделью
Рис. 4. Результат аппроксимации сигнала с шумом экспоненциальной моделью (метод Прони [1])
Для сравнения на рис. 4 представлены результаты аппроксимации этого сигнала модифицированным методом наименьших квадратов Прони [1], вычисляющим параметры модели (3) аналитически. Количество комплексных экспонент в модели Р = 21. Относительная ошибка аппроксимации (1) получилась равной
0.89, т.е. значения параметров, вычисленных методом Прони, соответствуют одному из локальных экстремумов приведенного в данной статье алгоритма.
Относительная ошибка аппроксимации тестовой последовательности Марпла [1], полученная программой при 10 поисковых параметрах, равна 0.0144 (рис. 5). Относительные ошибки частот не превышают 0.3%, относительные ошибки амплитуд - 3%. Время вычисления не превышает 10 минут (основная часть времени затрачивается на поиск более слабых компонент).
Данные исследования показали, что приведенный в статье алгоритм не чувствителен к аддитивному шуму и позволяет обнаружить в анализируемом сигнале как слабые частотные компоненты, так и близкие по частоте.
4. Частотно-временное распределение квазигармонического сигнала
Для анализа динамики частотных компонент квазигармонического сигнала вычислим оптимальные параметры синусоидальной (или экспоненциальной) модели приведенным алгоритмом в скользящем окне. Полученные значения амплитуд отображаются на плоскости время - частота. По оси абсцисс откладывается номер отсчета сигнала, соответствующий центру скользящего окна, по оси ординат - линейная цифровая частота, по оси аппликат на частотах гак/2тс отображаются значения амплитуд Ak оттенками серого цвета (белому цвету соответствует нулевая амплитуда, черному - максимальная амплитуда всего распределения).
Рассмотрим частотно-модулированный сигнал, цифровая частота которого меняется по синусоидальному закону от 0.01 до 0.02
Si = cos(2rc0.015+ 2.5cos(2ra'0.002))+
+ j sin (2^0.015+2.5cos(2ra'0.002)),
i = 0, 1, 2, ..., 999.
Рис. 5. Результат аппроксимации тестовой последовательности Марпла
Рис. 6. Частотно-временное распределение квазигармонического сигнала без шума
Рис. 7. Частотно-временное распределение квазигармонического сигнала с шумом
На рис. 6 изображено частотно-временное распределение данного сигнала, вычисленное с размером скользящего окна 50 отсчетов и двумя поисковыми параметрами. Количество итераций поиска оптимальных параметров в каждом положении окна равно 20.
На рисунке хорошо видна частотная компонента, закон изменения частоты которой полностью соответствует заданному (кроме начальных и конечных отсчетов, что связано с захватом окном краев сигнала).
При увеличении размера окна до 100 отсчетов появляется ошибка по частоте, которая не превышает 3%.
При добавлении к реальной и мнимой частям сигнала белого гауссова шума с дисперсией, равной 1, распределение, полученное с помощью скользящего окна в 50 отсчетов, становится размытым. Увеличение размера окна до 100 отсчетов позволяет выделить частотную
компоненту и качественно оценить закон изменения частоты (рис. 7).
Рассмотрим пример вычисления данным алгоритмом частотно-временного распределения реального сигнала со слабой частотной модуляцией на фоне мешающих факторов.
Анализируемый сигнал был получен на лабораторной установке «Маятник» [5], состоящей из маятника, в качестве груза которого используется металлическая пластина, хорошо отражающая звук, звуковых колонок, микрофона и персонального компьютера. Микрофон располагается между звуковыми колонками, излучающими синусоидальный сигнал частотой 4410 Гц, и маятником так, чтобы звук, излучаемый колонками, отражаясь от качающейся пластины, попадал в него. Сигнал, полученный с микрофона, оцифровывался звуковым контроллером компьютера с частотой дискретизации 44100 Гц.
Timej samples
Рис. 8. Осциллограмма квадратурного сигнала, полученного на установке «Маятник»
14.000
Timej samples
Рис. 9. Частотно-временное распределение сигнала, полученного на установке «Маятник»
Перед вычислением спектрально-временного распределения выполнялась предварительная обработка сигнала, которая заключается в фильтрации полосовым фильтром (с частотами среза 4 кГц и 5 кГц), квадратурном синхронном преобразовании и децимации (с шагом 50).
На рис. 8 изображены реальная (внизу) и мнимая (вверху) части сигнала, полученные после предварительной обработки.
На рис. 9 изображено спектрально-временное распределение, вычисленное с размером скользящего окна 100 отсчетов и одним поисковым параметром (еще 1 параметр с нулевой частотой и 1 с противоположной частотой добавляются в целевой функции). По оси абсцисс отложено время в отсчетах, по оси ординат -частота в герцах.
На спектрограмме хорошо видны две частотные компоненты: одна соответствует постоянной составляющей, а вторая - доплеров-ской частоте звуковой волны, отраженной от
движущейся пластины маятника. Доплеровская частота изменяется по синусоидальному закону (с амплитудой ~13 Гц или 0.015 в цифровой частоте) с периодом 750 отсчетов (0.85 с) и небольшим затуханием, что соответствует движению реального маятника.
Выводы
Описанный в данной статье алгоритм предназначен для оценки оптимальных параметров синусоидальной (2) и экспоненциальной (3) моделей сигнала, заданного на коротком отрезке.
Алгоритм не чувствителен к белому гауссо-вому шуму, что подтверждается примерами.
На примерах показана возможность использования данного алгоритма для вычисления частотно-временного распределения ква-зигармонических сигналов со слабой частотной модуляцией. Малый размер скользящего окна позволил увеличить разрешение по вре-
мени, при этом разрешение по частоте осталось высоким.
Недостатком алгоритма является длительное время оценки спектра.
Работа поддержана конкурсным контрактом КД НК-21П с Федеральным агентством образования и грантом № 228319 Европейского союза в рамках проекта EuroPlanet RI FP7.
Список литературы
1. Марпл-мл. С.Л. Цифровой спектральный анализ и его приложения. М.: Мир, 1990. 584 с.
2. Allen R.L., Mills D.W. Signal Analysis: Time, Frequency, Scale, and Structure. Wiley-IEEE Press, 2004. 966 p.
3. Cohen L., Time-frequency distributions - a review // Proceedings of the IEEE. 1989. V. 77. N. 7. P. 941-981.
4. Канаков В.А., Лупов С.Ю., Фрадкина Е.П. и др. Применение спектрально-временного анализа для исследования интерферометрических данных // Труды (девятой) Научной конференции по радиофизике «Факультет - ровесник Победы». 7 мая 2005 г. / Под ред. А.В. Якимова. Н. Новгород: ТАЛАМ, 2005. С. 117-118.
5. Лупов С.Ю., Фрадкина Е.П. Лабораторная установка для учебного курса «Цифровая обработка
сигналов» // Открытое образование. 2009. № 5. С. 30-34.
6. Дмитриев Е.В. Обнаружение и различение сигналов в их аддитивной смеси путем расчета и анализа ее естественного спектра // Физика волновых процессов и радиотехнические системы. 2008. Т. 11. № 2. С. 61-66.
7. Fradkina E.P. Short signals spectrals estimation method // Proceedings of the 12th Scientific Conference on Radiophysics dedicated to the 90th Anniversary of M.M. Kobrin’s Birth (Nizhny Novgorod, May 7, 2008) / Ed. by A.V. Yakimov, S.M. Grach. Nizhny Novgorod: TALAM Press, 2008. P. 272-73.
8. http://www.nag.co.uk/numeric/Fl/manual/pdf/E04 /e04ccf.pdf (дата обращения: 08.09.2010).
9. Лупов С.Ю., Фрадкина Е.П. Модификация функции LabVIEW «Downhill Simplex nD» и тестирование её на примере аппроксимации тестовых сигналов суммой синусоид // Труды Международной научно-практической конференции «Образовательные, научные и инженерные приложения в среде LabVIEW и технологии National Instruments». М.: Изд-во РУДН, 2006. C. 278-281.
10. http ://www. nag.co. uk/numeric/fl/manual 19/pdf/E 04/e04jaf_fl18.pdf (дата обращения: 08.09.2010).
11.Nelder J.A., Mead R. A simplex method for function minimization // Computer Journal. 1965. V. 7. P. 308-313.
12. Гилл Ф., Мюррей У., Райт М. Практическая оптимизация: Пер. с англ. М.: Мир, 1985. 509 с.
EVALUATION OF OPTIMAL PARAMETERS OF DISCRETE SIGNAL SEGMENT EXPONENTIAL AND SINUSOIDAL MODELS
S.Yu. Lupov, A.M. Serebryakov, E.P. Fradkina
An algorithm has been proposed to evaluate optimal parameters of sinusoidal and exponential models of a discrete signal segment. It has been shown that the algorithm can be applied to evaluate the time-frequency distribution of a quasi-harmonic signal. An example is presented of time-frequency distribution calculation of a real signal obtained on the laboratory facility “Mayatnik” (pendulum).
Keywords: spectrum estimation, time-frequency distribution, optimization.
