Оценка несущей частоты случайной последовательности импульсов методом Прони Текст научной статьи по специальности «Электротехника, электронная техника, информационные технологии»

УДК 519.676

ОЦЕНКА НЕСУЩЕЙ ЧАСТОТЫ СЛУЧАЙНОЙ ПОСЛЕДОВАТЕЛЬНОСТИ ИМПУЛЬСОВ МЕТОДОМ ПРОНИ

А. Л. Шестаков, А. С. Семенов, О. Л. Ибряева

CARRIER FREQUENCY ESTIMATION FOR RANDOM PULSE TRAIN USING PRONY'S METHOD

A. L. Shestakov, A. S. Semenov, O. L. Ibryaeva

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

Ключевые слова: сегментированный метод Прони, оценка частоты, затухающая синусоида

This paper presents an analysis by means of mathematical model approach for system response to the random pulse train using the segmented Prony's method. The system response is the damped sinusoids sequence, and sinusoids can be partially overlapped. The observational interval is segmented in short segments, and Prony's method is applied to each of them. To characterize the method performance the carrier frequency variance estimate is used.

Keywords: segmented Prony's method, frequency estimation, damped sinusoid

Введение

Задача оценки параметров комплексного экспоненциального сигнала является одной из самых распространенных в различных областях техники. Это связано с тем, что отклик линейной системы на импульсное воздействие является суммой именно таких сигналов, т.е. оценив параметры сигналов на выходе системы, можно решить задачу идентификации системы и ее состояния. Попытка использования в этих целях преобразования Фурье дает не очень хорошие результаты, так как оно предназначено для оценки спектра сигнала, а не частоты и, кроме того, в классическом варианте не является статистически устойчивым. Метод Прони изначально был предназначен для оценки параметров линейной комбинации экспоненциальных функций [1], что позволяет сделать предположение о возможности его применения для диагностики линейных систем. Одним из допущений при разработке алгоритма обработки сигналов по методу Прони является предположение о точно известном времени начала сигнала, в связи с чем и при анализе этого метода делается такое же допущение [2]. В то же время, в случае действия на входе линейной системы случайного потока импульсных сигналов типа ¿-функции, это предположение нельзя считать верным. Таким образом, возникает задача модификации метода Прони для потока импульсов в виде комбинации комплексных экспонент, возникающих в случайные моменты времени, и оценки его характеристик для этих условий.

1. Модель сигнала

Входным возбуждающим воздействием в рассматриваемой модели является последовательность ¿-импульсов, возникающих в случайные моменты времени Г\,Г2,... ,т8,..., тм• Под действием каждого из ¿-импульсов в системе возникает затухающее колебание с частотой и = 2п/ и коэффициентом затухания а. Для простоты в данной работе мы ограничиваемся случаем одной собственной частоты датчика, т.е. считаем, что система представляет собой резонансный контур с собственной частотой /.

Амплитуда колебания А3 = А (т8) является нормальной случайной величиной с математическим ожиданием равным нулю и среднеквадратичным отклонением а. Величина А3 и есть «высота> ¿-импульса, т.е. входное воздействие и(г) имеет вид: и(г) = ^М=1 А3 5 (г — т8). Выходной сигнал представляет собой последовательность затухающих синусоид, возникающих в моменты времени т8:

м

Б(1) = £ А5 • V(г — т5) е

(а+г 2п/)

(1)

8=1

где V (г) =

1, г > о,

Фаза р (т8) - равномерная случайная величина на — ", 2].

0, г < о. -—гу», — — --------ь 2 > ^

Пример входного потока импульсов и соответствующего ему выходного сигнала приведен на рис. 1. (В этом и последующих примерах статьи мы брали частоту дискретизации равной 250 кГц, собственную частоту / = 32 кГц, коэффициент затухания а = —20 кГц и а = 0, 5.) Заметим, что синусоиды могут «наползать> друг на друга.

Рис. 1. Входное воздействие и выходной сигнал

В общем случае, к выходному сигналу должен быть добавлен шум, который считаем белым, нормально распределенным N (0, а\). Отношение сигнал/шум (SNR) для данной модели сигнала будем вычислять по формуле:

БЖ = 10 ^

1 М 1 ие+М 2

М ^ N , ^ у2

8=1 к=иа

а

Здесь М - число 5-импульсов на рассматриваемом временном отрезке сигнала, N - число отсчетов, за которое синусоида затухает до уровня 5 % от энергии в момент возникновения импульса, п8 - номер отсчета, соответствующий времени т8, у к - значения выходного сигнала в отсутствии шума.

Пример зашумленного сигнала для SNR=10 дБ приведен на рис. 2. Там же в увеличенном масштабе показан вид отклика на один импульс.

О 0.002 0.88 Щ

Игле, гее (¡те, 10"4

Рис. 2. Зашумленный сигнал с 8МИ=10 дБ

2. Предварительные замечания

Метод Прони обычно используется для сигналов, подобных приведенному на рис. 2. С точки зрения теории оценки параметров это означает, что задача обнаружения решена, и начало сигнала известно точно, т.е. требуется оценить только амплитуду А сигнала, постоянную затухания а, частоту и и фазу ф. В случае нескольких синусоид, помимо параметров каждой, дополнительно требуется оценить и их число. Для практики обычно наибольший интерес представляют частота и постоянная затухания. Именно в такой постановке исследуется метод Прони в [2, 3].

Желание использовать метод Прони для представленного на рис. 2 сигнала (1) с шумом, приводит нас к следующим вопросам:

• какой длительности необходимо выбирать окно обработки,

• как устанавливать окно обработки (как определить начало импульса),

• что происходит, если начало окна обработки не совпадает с началом импульса,

• как работает алгоритм, если импульсы перекрываются.

Таким образом, при сигнале вида (1) мы сталкиваемся не только с задачей оценки параметров, но и с задачами обнаружения и различения сигналов.

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

Такой способ обработки носит название сегментированного метода Прони и используется, например, в [4]. Впервые же термин «сегментированный метод Прони> был, насколько нам известно, применен в [5]. Длина сегментов в этой работе не была, вообще говоря, одинаковой, кроме того, сдвиг осуществлялся на длину целого сегмента, с переходом к следующему. Важно отметить, что модельный сигнал (1) отличается от нестационарных сигналов,

рассматриваемых в работах [4, 5]. В указанных работах нестационарность определяется изменением частоты на интервале наблюдения. Напротив, в рассматриваемом случае частота на исследуемом отрезке постоянна, но сам полезный сигнал присутствует не всегда.

3. Применение сегментированного метода Прони для оценки параметров потока экспоненциально затухающих синусоид

Как уже отмечалось, мы будем использовать модификацию метода Прони, состоящую в том, чтобы обрабатывать последовательно отрезки сигнала длины N << Ь (Ь = 100 000 -число отсчетов всего сигнала), сдвигаясь на единицу. Для анализа составляющих сегментов будем использовать совместно сингулярные методы Прони вперед и назад [1]. Выделим несколько основных этапов.

Сначала по N отсчетам сигнала определяются параметры линейного предсказания вперед и назад. Затем из них формируются характеристические полиномы вперед А (г) и назад В (г), обладающие следующим свойством: корни полиномов А (г), В (г) есть соответственно:

$ = е(ак +г2п!к )Т гк = е(-аk+i2пfk )Т (2)

(Здесь 1 < к < р, р - число комплексных экспонент, присутствующих в сигнале, ак, /к - коэффициент затухания и частота к-й экспоненты, Т - период дискретизации сигнала.) Наконец, зная корни полиномов, можно найти частоты /к.

Хорошо известно (см., например, [1]), что улучшению результатов способствует завышение числа истинных экспонент, присутствующих в сигнале. В этом случае, при анализе нулей обоих полиномов А (г) и В (г), возможно их разделение на истинные нули сигнала и нули шума. Мы также будем использовать завышенный порядок предсказания для обработки нашего сигнала.

Заметим далее, что, поскольку одна затухающая синусоида с частотой / соответствует двум комплексным экспонентам с частотами / и —/, то, ясно, что в этом случае 1/г^ будут являться корнями многочлена В (г) и, наоборот, 1/г| будут корнями А (г). (Конечно, такое равенство корней в данном случае из-за наличия шума будет справедливо лишь приближенно.) Тогда в дальнейшем можно брать, например, корни многочлена В (г), лежащие в верхней полуплоскости и по модулю большие единицы и корни многочлена А (г) из нижней полуплоскости, по модулю меньшие единицы. Среди всех этих корней мы будем искать обратные друг к другу.

Здесь важно отметить серьезную трудность, которая возникает при разработке программного алгоритма: из-за наличия шума нам не удастся найти такие корни, что

г\,--а = 0, и данное равенство будет справедливо лишь приближенно. Будем поэтому

к Х1

считать корни г\,, г" обратными друг к другу, и не удовлетворяющими точному соотношению выше только по причине шума, если

Л 1

гк г1

< £, (3)

где е > 0 - параметр, значение которого мы задаем сами.

На рис. 3 построены корни многочлена В (г) (из верхней полуплоскости, большие по модулю единицы) и корни многочлена А (г) (из нижней полуплоскости, меньшие по модулю единицы), которые были получены при обработке 100 000 отсчетов модельного сигнала с SNR=10 дБ сегментированным методом Прони. Для обработки каждого сегмента длиной N = 50 отсчетов был использован завышенный порядок предсказания р =12 (вместо истин-

ного порядка, равного двум). На правой половине рис. 3 показаны только те корни, которые удовлетворяют неравенству (3) при е = 0, 01. Они группируются возле положений истинных нулей сигнала. Именно по таким корням, используя соотношения (2), мы находим частоту сигнала. Учитывая, что на основе неравенства (3) определяется наличие сигнала, можно назвать параметр е радиусом зоны обнаружения (РЗО) сигнала.

-1 -0.5 0 GÏS 't -1 D G.I '1

Рис. 3. Корни полиномов линейного предсказания

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

Рис. 4. Зависимость среднеквадратического отклонения от РЗО

Ясно, что с одной стороны, меньшее значение е дает более точную оценку частоты. Однако, выбирая малые значения е, мы возможно не получим никакой (даже «плохой>) оценки частоты, так как не найдется ни одной пары столь близких корней. Кроме того, при малых значениях е и плохом соотношении сигнал/шум, может происходить случайное попадание отдельных измерений в РЗО, что приведет к аномальным выбросам оценок. Поэтому е не может быть выбрано слишком малым.

Исходя из приведенных соображений, примем далее е=0,01, т.к. при этом значении заканчиваются осцилляции СКО на рис. 4 при SNR=1 дБ, которые могут явиться признаком наличия аномальных оценок. При данном е обработаем сегментированным методом Про-ни 100 000 отсчетов нашего модельного сигнала, сдвигаясь каждый раз на один отсчет и применяя сингулярные методы Прони вперед и назад с порядком предсказания р =12 для обработки составляющих сегментов длиной N = 50 отсчетов. В результате мы обнаружим целый набор найденных частот сигнала. (Разумеется, частоты найдутся не на всех сегментах.)

С помощью полученного набора частот построим гистограмму на рис. 5. Здесь по оси абсцисс отложены значения частот от 0 до 125 кГц (половина частоты дискретизации), а по оси ординат число появлений данной частоты. Таким образом, как и в случае спектра Фурье, этот рисунок показывает вклад каждой из частот в обработанный сигнал. Мы будем называть поэтому такую гистограмму спектром Прони. Заметим, что в [1] принято другое определение спектра Прони. Оно не подходит нам, поскольку рассматриваемый сигнал не укладывается в рамки допущений относительно вида колебаний, сделанных в [1]. По этой же причине в статье [5] для сигнала, частота которого может резко меняться со временем, используется иной способ вычисления спектра Прони. Нам в этой статье будет удобно вычислять спектр Прони, как описано выше.

Рис. 5. Спектр Прони

На рис. 6 построена зависимость среднеквадратического отклонения а оценки частоты от соотношения сигнал/шум.

0|-1-1-1-1-1-

О 5 10 15 20 25 30

БШ. №

Рис. 6. Зависимость среднеквадратического отклонения а оценки частоты от соотношения сигнал/шум

На рис. 7 приведена зависимость отклонения 5 = |/ — /| мат. ожидания частоты / (от истинного значения f = 32 000 Гц) от соотношения сигнал/шум.

140-,-.-.-,-.-

х

СО

тз

0 5 10 15 20 25 30 ЭЫ«, ЙВ

Рис. 7. Зависимость отклонения мат. ожидания частоты (от истинного значения) от соотношения сигнал/шум

4. Применение преобразования Фурье для оценки параметров потока экспоненциально затухающих синусоид

Для анализа перспектив использования метода Прони в задаче оценки параметров экспоненциального сигнала, рассмотрим характеристики преобразования Фурье при обработке сигнала (1), описанного выше. Чтобы повысить достоверность сравнения методов Прони и Фурье, необходимо максимально совместить условия обработки.

С этой целью метод Фурье реализуется в виде модификации Уэлча с окном обработки N = 50 отсчетов и смещением окна обработки на один отсчет, что аналогично параметрам сегментированного метода Прони, рассматриваемого ранее.

Для оценки влияния размеров окна на характеристики преобразования Фурье все расчеты были проведены также для размеров окна N = 1000. На рис. 8 показан спектр чистого сигнала и сигнала с SNR=1 дБ для значений N = 50, 1000.

Рис. 8. Спектр чистого модельного сигнала и сигнала с 8МИ=1 дБ, полученный методом Уэлча для окна размеров N = 50, 1000

После получения спектра на каждом окне обработки найдем частоту, имеющую наибольшую мощность и повторим эту операцию при каждом смещении окна. Тогда можно построить гистограмму, аналогичную приведенным на рис. 5 и показанную на рис. 9 N = 50) и рис. 10 N = 1000).

without noise

2.5 3 3.5

frequency, Hz

2.5 3 3.5

frequency, Hz

Рис. 9. Гистограмма частот для случая N = 50

Эти гистограммы позволяют вычислить необходимые статистические характеристики метода. Результаты расчета при различных значениях шума приведены на рис. 11 (зависимость среднеквадратического отклонения от SNR) и рис. 12 (зависимость отклонения от мат. ожидания от SNR).

without noise SNR--1 dB

frequency, Hz 4 frequency, Hz

Рис. 10. Гистограмма частот для случая N = 1000

SNR, dB

Рис. 11. Зависимость среднеквадратического отклонения частоты от соотношения сигнал/шум

2500 -N=50 --- N=1000

2000 —-—

1500

iooo

.S00

0

О -'-т----'----'-----

□ 5 10 15 20 Ж 30

И

Рис. 12. Зависимость отклонения мат. ожидания частоты (от истинного значения) от соотношения сигнал/шум

Заключение

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

1. Модифицированный метод Прони позволяет достигнуть погрешности оценки частоты 1 % уже при соотношении сигнал/шум около 15 дБ и 50 отсчетов окна наблюдения. Заметим, что преобразование Фурье не позволяет добиться такого результата даже при соотношении сигнал/шум 30 дБ и 1000 отсчетов окна наблюдения. Таким образом, метод Прони имеет большие перспективы при решении задачи оценки частоты экспоненциально затухающего сигнала и, в частности, может быть использован для оценки состояния линейных систем.

2. В данной статье впервые введено понятие радиуса зоны обнаружения сигнала. Для практического использования можно рекомендовать значение параметра РЗО около 0,01 независимо от соотношения сигнал/шум и иных параметров сигнала.

3. В целях разработки рекомендаций для решения практических задач в дальнейших работах необходимо рассмотреть ряд вопросов, имеющих существенное влияние на характеристики метода, а именно: 1) оценить возможность работы метода при многочастотных импульсах; 2) рассмотреть влияние перекрытия импульсов на характеристики метода; 3) исследовать влияние отношения длительности окна к длительности импульса на точность метода; 4) рассмотреть параметры оценки постоянной затухания; 5) оценить возможность обработки массива отсчетов, полученных при сегментарной обработке, для повышения точности метода.

Работа выполнялась в рамках госконтракта НК-166П на поисковую научно-исследовательскую работу «Разработка и исследование методов оценки состояния преобразователя давления в ходе технологического процессам.

Литература

1. Марпл-мл., С. Л. Цифровой спектральный анализ и его приложения: пер. с англ. / С. Л. Марпл-мл. - М.: Мир, 1990. - 584 с.

2. Кузнецов, Ю. В. Сравнительная характеристика алгоритмов оценки параметров резонансной модели объектов / Ю. В. Кузнецов, В. Ю. Щекатуров, А. Б. Баев // Вестн. МАИ. - 1997. - Т. 4, №2. - С. 70 - 76.

3. Кузнецов, Ю. В. Использование метода Прони и его модификаций при оценке параметров резонансной модели / Ю. В. Кузнецов, А. Б. Баев // Будущее авиации и космонавтики: сб. тез. ст. науч.-исслед. работ студентов. - М.: Изд-во МАИ, 1998. -С. 46 - 49.

4. Кухаренко, Б. Г. Технология спектрального анализа на основе быстрого преобразования Прони / Б. Г. Кухаренко // Информ. технологии. - 2008. - №4. - С. 38 - 42.

5. Barone, P. The segmented Prony method for the analysis of non-stationary time series / P. Barone, E. Massaro, A. Polichetti // Astron. Astrophys. - 1989. - V. 209. - P. 435 - 444.

Кафедра информационно-измерительной техники, Южно-Уральский государственный университет [email protected]

Поступила в редакцию 10 сентября 2009 г.