Оценка частоты гармонических сигналов на основе анализа амплитудно-частотной характеристики процесса

Размещено на http://www.allbest.ru/

Оценка частоты гармонических сигналов на основе анализа амплитудно-частотной характеристики процесса

Бордюков Антон Геннадьевич, аспирант Санкт-Петербургского Государственного Электротехнического Университета им. В.И. Ульянова (Ленина), инженер 2-й категории ЗАО «Моторола ЗАО»

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

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

Рассмотрим дискретизированную реализацию xn изучаемого процесса x, состоящую из N последовательных отсчетов x(nt), где n = 0,…, N-1; t - временной интервал, через который получены точки дискретизации. Процесс x представляет собой сумму белого шума и нескольких гармонических компонент, имеющих постоянную амплитуду и частоту в пределах времени измерения Nt. Задачей разрабатываемого алгоритма является определение частоты рассматриваемых гармонических компонентов.

В основе получения АЧХ дискретного сигнала лежит дискретное преобразование Фурье (ДПФ) X(fk), вычисляемое по формуле [1]:

,(1)

где wn - весовая функция окна, применяемая к входным данным; kw - нормировочный коэффициент, компенсирующий искажения, вызванные функцией окна. Коэффициент kw может быть определен следующим образом [1]:

.

При этом частота fk принимает дискретные значения, определяемые по формуле:

, k=0, 1, …, N-1

АЧХ ДПФ выражается формулой:

,(2)

где XRe(fk) и XIm(fk) - вещественная и мнимая часть ряда X(fk), получаемого в результате выполнения ДПФ.

Характеристики ДПФ, оказывающие наибольшее влияние на точность определения частоты гармонических сигналов.

Основным параметром, определяющим точность оценки частоты, является шаг дискретизации fs. Для дискретизированного процесса с частотой дискретизации fs шаг fs = fs/N. Данный параметр удобно использовать, как единицу измерения абсолютной погрешности.

ДПФ конечной последовательности данных может быть представлено как перемножение процесса, заданного на бесконечном интервале времени, с прямоугольным «временным окном» wn (взвешивающей, или весовой, функцией), которое можно определить следующим образом [2]:

Перемножению во временной области соответствует свертка в частотной области, то есть , где Wk и Xk - ДПФ оконной функции и исследуемого процесса соответственно, а символом обозначена операция свертки [2]. Весовую функцию Wk принято называть «спектральным окном» [2]. Формулу для прямоугольного спектрального окна можно выразить следующим образом [3]:

,

где величина k определена в любой точке от 0 до , причем изменение k на 1 соответствует изменению частоты на один шаг дискретизации fs.

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



График АЧХ свертки прямоугольного и треугольного спектрального окна с ДПФ гармонического сигнала частотой 6.5fs, полученный по формулам (1, 2), приведен на рисунке 1.

Рис. 1 - Иллюстрация параметров, связанных с эффектом просачивания

Кратко рассмотрим основные параметры спектральных окон, оказывающих влияние на результирующую АЧХ сигнала. Данные параметры проиллюстрированы на примере треугольного окна (Рис. 1).

1. Гребешковое искажение SL (scalloping loss) (дБ) определяет, насколько уменьшится амплитуда сигнала относительно истинного значения амплитуды в точках дискретизации в худшем случае, то есть когда гармоника попадает точно между двумя точками дискретизации.

2. Потери при обработке PL (processing loss) (дБ). Не прямоугольные окна за счет уменьшения амплитуды сигнала на краях последовательности имеют меньшую амплитуду основного лепестка, нежели исходный сигнал, что в конечном итоге ведет к снижению фактического значения амплитуды в точках дискретизации. Данное снижение амплитуды называется потерями при обработке.

3. Потери в худшем случае (worst case process loss) (дБ) складываются из потерь при обработке PL и гребешкового искажения SL и, по сути, являются снижением отношения сигнал/шум, вызванного взвешиванием и наихудшим расположением гармонической составляющей.

4. Максимальный уровень первого бокового лепестка HSLL (highest side lobe level) (дБ) - отношение амплитуды первого бокового лепестка к амплитуде главного лепестка.

5. Асимптотическая скорость спада боковых лепестков SLFO (side lobe fall off) (дБ/октава) определяет скорость спада амплитуды боковых лепестков от лепестка к лепестку.

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

Ряды для случаев прямоугольного и треугольного окна приведены в таблице 1. В первой строке указано значение k, а ячейки второй строки показывают отношение истинной амплитуды гармоники (0 дБ) к амплитуде в точке k.

Таблица 1 - Некоторые значения АЧХ функций прямоугольного и треугольного спектрального окна

Для прямоугольного окна
k, fs	0.5	1.43	2. 46	3.47	4.48	5.48	6.49	7.49	8.49	9.49
Wl(k), дБ	-3.9	-13.3	-17.8	-20.8	-23	-24.7	-26.2	-27.4	-28.5	-29.5
Для треугольного окна
k, fs	0.5	1	2.86	-	4.92	-	6.94	-	8.96	10.96
Wl(k), дБ	-3.07	-9.1	-27.8	-	-36.9	-	-42.8	-	-47.2	-50.7

 

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

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

, (3)

где A12 - отношение амплитуд гармоник, A12 = A1/A2; f1 и f2 - истинное значение частот гармоник.

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

Алгоритм обнаружения гармонических составляющих.

Алгоритм обнаружения построен на анализе относительных амплитуд АЧХ ДПФ, выраженных в децибелах, где значению АЧХ, равному 1, соответствует уровень 0 децибел.

Для принятия решения об обнаружении гармонической составляющей по АЧХ необходимо ввести соответствующий критерий принятия решения. В основу данного критерия положим принцип превышения амплитудой порогового значения X0. Кроме того, для улучшения характеристик обнаружения целесообразно учитывать особенности воздействия спектральных окон. Для этого применим дополнительный критерий, суть которого заключается в превышении значением амплитуды в точке дискретизации значений в соседних точках дискретизации на величину перехода X децибел. Децибелы удобно использовать по причине изменения значений в широком диапазоне амплитуд, а также в силу относительной природы спектральных окон. С учетом описанных принципов был получен алгоритм обнаружения гармонических составляющих процесса. Упрощенно данный алгоритм можно выразить следующей последовательностью шагов.

1. Предварительно проверяем все значения АЧХ: Если значение в точке i меньше значения X0 (децибел), то будем считать, что в данной точке не содержится значимых гармонических составляющих и будем принимать амплитуду в этой точке равной X0.

2. Последовательно рассматриваем пару соседних точек дискретизации по частоте i и i+1. Если значение амплитуды в точке i+1 превышает значение в точке i более чем на величину X децибел, то запоминаем точку i+1, так как в точке fmL=i+1 потенциально может быть зафиксирована гармоническая составляющая. Повторяем шаг 2 до тех пор, пока значение в точке i+1 больше или равно значению в точке i. В случае, если на последующих итерациях амплитуды в соседних точках i и i+1 отличаются более чем на значение X, то фиксируем точку fmL=i+1, считая её наиболее близкой к истинному положению гармоники. В случае если значение в точке i+1 меньше значения в точке i и значение fmL зафиксировано, переходим к шагу 3 и начинаем поиск необходимого спада АЧХ.

3. Последовательно рассмотрим пару соседних точек дискретизации по частоте i и i+1. Если значение в точке i превышает значение в точке i+1 на величину X децибел, то запоминаем точку fmH=i. Если значение в точке i меньше значения в точке i+1, тогда возвращаемся к шагу 2. Если точка fmH зафиксирована, то определяем положение гармоники как fm = (fmL+fmH)/2 и переходим к поиску следующей гармоники, увеличивая величину m на 1 и возвращаясь ко второму шагу алгоритма.

Данный алгоритм позволяет определить частоту отдельной гармоники при условии отсутствия помех и отсутствия влияния соседних гармоник с абсолютной погрешностью f=1/2fs. Однако в условиях реальных измерений данное значение может существенно изменяться в зависимости от параметров входного сигнала и выбранных параметров алгоритма.

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

Важнейшей характеристикой, оказывающим влияние на оба параметра критерия, является минимальная амплитуда сигнала A, выраженная в децибелах относительно истинной амплитуды гармоники, позволяющая зафиксировать гармоническую составляющую при любом её расположении. Данная величина является функцией X. Математически строгий расчет данной величины является сложной задачей, однако аналитически минимальную амплитуду, необходимую для обнаружения гармоники, можно оценить с помощью формулы:

(4)

Выражение (4) применимо для следующих значений величины перехода X:

,(5)

Если значение X-(PL+WI(1)) и при этом выполняется условие (4,а), то для оценки параметра A используется следующая формула:

(6)

Большое значение минимальная амплитуда сигнала имеет для определения величины порога X0.Также на выбор порога X0 оказывает влияние диапазон амплитуд гармоник, которые требуется обнаружить в ходе анализа АЧХ, с одной стороны, и свойства используемой функции окна с другой стороны. Так, даже в простейшем для обнаружения случае, когда сигнал представляет собой единственную гармонику, и помехи отсутствуют, применение большинства частотных окон дает результирующую амплитудно-частотную характеристику, содержащую ложные пики. Максимальная амплитуда ложных пиков может достигать уровня первого бокового лепестка. Учитывая эффекты потерь при обработке и в случае худшего расположения гармоники между точками дискретизации, требования к выбору параметра X0 можно выразить следующей формулой:

,(7)

где Amax, Amin (дБ) - максимальная и минимальная амплитуда АЧХ, которая может быть обнаружена.

Подставим (5) в (7) и выразим отношение амплитуд Amax-Amin:

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

Небольшое число спектральных окон (прямоугольное окно, окно Ханна с =2, окно Кайзера, окно Блэкмэна, окно Тьюки с >1) позволяет избежать образования ложных пиков в описанном случае. Это связано с тем, что данные спектральные окна имеют форму АЧХ монотонно убывающую по мере удаления точек дискретизации от истинного положения гармоники при любом расположении гармоники в пределах одного интервала дискретизации. Это объясняется тем, что ширина боковых лепестков равна одному интервалу дискретизации по частоте. Также большое значение имеет формой боковых лепестков. Для таких спектральных окон справедливо следующее соотношение:

(8)

В дальнейшем будем называть такие спектральные окна монотонными.

Формулы (7), (8) показывают, что значение порога X0 зависит от величины перехода X. Увеличение величины перехода приводит к понижению максимального допустимого порога X0, что в свою очередь приводит к повышению вероятности ложной тревоги. Следовательно, большую важность приобретает минимально допустимый уровень величины перехода X.

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

,

где Anoise - максимальная амплитуда шума с вероятностью P, выраженная в децибелах.

Минимальный уровень перехода X также определяется на основании информации о шуме, имеющемся в изучаемом процессе. Так, для минимизации вероятности ложной тревоги необходимо рассмотреть худший случай, когда боковой лепесток имеет амплитуду равную пороговому значению X0, а аддитивный шум усиливает данную величину при неизменном значении соседних точек дискретизации. Используя формулу (8) для определения минимального значения X0 получаем следующее отношение:

где amin - абсолютное значение минимальной амплитуды АЧХ, которую необходимо обнаружить в спектре, эквивалентное относительной амплитуде Amin; anoise - абсолютное значение, которое с вероятностью P не будет превышено амплитудой шума.

Подставив значение параметра A(X) из формулы (4) и выразив параметр X, получаем следующий результат:

,

где (9)

Левая граница диапазона устанавливает минимальное значение для величины перехода X, при котором отсутствие ложной тревоги по причине воздействия белого шума будет гарантировано.

Не сложно показать, что при условии большого количества равномерно распределенных гармоник критерий не позволит обнаружить отдельные гармонические составляющие с приемлемой точностью. Придельный случай, когда существует возможность обнаружить 2 гармонические составляющие с помощью описанного критерия, определяется разностью частот гармоник f12=f1-f2. Данная величина является важным параметром, так как от неё напрямую зависит вероятность обнаружения разных по амплитуде гармоник. Для определения предельного значения разности частот между гармониками необходимо обладать априорными данными о диапазоне амплитуд гармонических составляющих A=Amax-Amin, которые требуется обнаружить с помощью описанного критерия, а также выбрать значение X на основании описанных выше рекомендаций. Кроме того, большое значение для расчета погрешности имеет форма используемого спектрального окна.

Учитывая наихудшее взаимодействие фаз (5), f12 можно оценить по следующей формуле:

,(10)

где в качестве величины kx предполагается ближайший локальный максимум. Значение kx определяет расстояние между гармониками f12, выраженное в интервалах дискретизации fs.

Выражение (10) дает несколько завышенную оценку, так как при выводе формулы принято допущение, что WI(kx-1)=WI(kx). Справедливость данного допущения зависит от полученной величины kx, а также выбранного спектрального окна. Так, для случаев окон Бартлетта, Бартлетта-Хэннинга, Бохмана [4] или треугольного окна при kx 3 можно считать WI(kx-1)=WI(kx). в силу того, что для данных окон характерны боковые лепестки шириной 2 интервала дискретизации. Для случая окна Чебышева все боковые лепестки имеют одинаковую максимальную амплитуду, что позволяет считать выражение (10) точным для любого kx большего ширины главного лепестка спектрального окна.

Сравнение спектральных окон.

Сравним результаты, получаемые с помощью выбранного критерия, для различных спектральных окон. Сведем в таблицу основные характеристики нескольких спектральных окон. Кроме того, приведем максимальные значения X, рассчитанное на основании условия (5), и значения f12 для X = 3дБ и для разного соотношения амплитуд A.

Таблица 2 - Характеристики спектральных окон

Окно	HSLL (дБ)	SLFO (дБ/окт)	kw (дБ)	SL (дБ)	PL (дБ)	WCPL (дБ)	X (дБ)	f12 (fs) при X=3, A=
								0 дБ	20 дБ	40 дБ
Прямоугольное*	-13	-6	0	3.92	0	3.92	6.5	2.5	27.5	>200
Треугольное	-27	-12	-4.77	1.82	1.25	3.07	4	2.86	6.95	-
Ханна, =2*	-32	-18	-4.28	1.42	1.76	3.18	2.9	2.5	3.5	7.5
Хэмминга	-43	-6	-4.02	1.78	1.32	3.1	3	2.5	3.5	-
Чебышева, 50дБ	-50	0	-4.16	1.7	1.42	3.12	2.7	2.3	2.3	-
Чебышева, 70дБ	-70	0	-4.95	1.25	2.1	3.35	2.4	3	3	3
Блэкмэна*	-58	-18	-5.18	1.1	2.37	3.47	2.2	3.5	3.5	4.5

Символом * помечены монотонные функции окна.

Не монотонные функции требуют выполнения условия (7), в результате чего невозможно определить значение f12 для наибольшего диапазона амплитуд для ряда таких окон.

Обращают на себя внимание крайне низкие значения X для большого числа спектральных окон, что делает для них практически непригодной формулу (4). В этом случае для расчета X необходимо воспользоваться формулой (5), что ведет к увеличению значения A(X).

Существенное увеличение значения A(X) ведет к ухудшению параметров f12 и X0, что может служить причиной увеличения количества ложных срабатываний при значительной амплитуде шума. Формула (5) была использована в расчетах величины f12 в таблице 2 для спектральных окон, значение X для которых меньше 3 дБ.

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

Имитационное моделирование.

Моделирование производилось в среде Matlab. В качестве входной реализации процесса рассматривалась последовательность точек длиной N=256. Данная последовательность содержала 2 синусоиды, имеющие произвольную начальную фазу, частоту и амплитуду. Значение амплитуды и частоты было постоянно в пределах одной реализации. Значение начальной фазы для каждой гармоники распределялось по равномерному закону в диапазоне от 0 до 2. Равномерное распределение также использовалось для частоты и амплитуды, причем диапазон амплитуд варьировался в зависимости от проводимого опыта. Так рассматривались случаи равных амплитуд, диапазона амплитуд от 0 до 20 дБ и от 0 до 40 дБ. Диапазон частот был ограничен таким образом, чтобы гармонические составляющие спектра находились на расстоянии не менее 2 интервалов дискретизации от частот 0 и fs/2.

На основании формул (1) и (2) с помощью одного из весовых окон, упомянутых в таблице 2, была получена АЧХ процесса. Затем значения амплитуд АЧХ приводились к децибелам и обрабатывались описанным алгоритмом. Результатом работы алгоритма являлись значения частоты, которые сравнивались с частотами гармоник, поступающими на вход функции. Приведенная последовательность повторялась 10 000 раз, после чего, были получены среднее и максимальное значение погрешности. Результаты моделирования сведены в таблицу 3.

сигнал фурье дискретный преобразование

Таблица 3 - Результаты моделирования

Окно	f12 (fs) при X=3, A=	f (fs) при X=3, A=	f (fs) при X=1, A=
	0 дБ	20 дБ	40 дБ	0 дБ	20 дБ	40 дБ	0 дБ	20 дБ	40 дБ
Прямоугольное*	2.6	13.7	117	0.19	0.21	0.81	0.23	0.24	0.39
Треугольное	2.67	3.07	122	0.14	0.16	3.9	0.21	0.21	3.9
Ханна, =2*	2.68	3.35	4.8	0.14	0.15	0.17	0.19	0.2	0.21
Хэмминга	2.7	3.15	55	0.14	0.15	0.23	0.2	0.21	0.24
Чебышева, 50дБ	2.67	3.2	39	0.14	0.15	0.17	0.19	0.2	0.23
Чебышева, 70дБ	3.2	3.4	4.2	0.15	0.16	0.17	0.18	0.19	0.2
Блэкмэна*	3.4	4.2	4.1	0.16	0.17	0.18	0.18	0.18	0.2

 

В качестве значения X0 применялась величина, отвечающая условию (7) и (8) в зависимости от типа применяемого окна.

Полученные результаты в целом соответствуют расчетным. Так значение f12 для отношения амплитуд A=40дБ для не монотонного спектрального окна, подтверждает формулу (7), а значение абсолютной погрешности f не превышает 0.5fs по модулю. Как аналитический расчет, так и имитационное моделирование, подтвердили, что наилучшие результаты в широком диапазоне амплитуд получаются с помощью окна Ханна.

Единственное несоответствие получено для величины f12 для A=0. Данная величина превышает расчетное значение для любой функции окна. Несоответствие можно объяснить тем, что при малых значениях f12 взаимодействию между собой подвергаются главные лепестки, в то время как вывод соотношения (10) осуществлялся при предположении взаимодействия боковых лепестков. Данное предположение косвенно подтверждается тем, что для окна Чебышева завышенными получаются значения для любого диапазона амплитуд. Так, для данного типа окна параметр f12 зависит только от размера главного лепестка, в то время как боковые лепестки не способны оказывать воздействие на предельное значение погрешности.

С целью проверки соотношения (9) к входному сигналу был дополнительно добавлен белый шум, дисперсия которого была выбрана таким образом, чтобы в результате применения формулы получался вещественный интервал значений для параметра D. После этого значение X, используемое при вычислениях, изменялось, чтобы в одном случае входить в указанный диапазон, а в другом выходить за пределы диапазона. Полученные результаты свидетельствовали о завышенной жесткости соотношения (8). Так резкое изменение погрешности при переходе через границу Xmin наблюдалось в случае использования в формуле (9) величины PL вместо WCPL. Данную корректировку можно объяснить тем, что условие (4), используемое для расчета формулы (9), соответствует наихудшему, реально не достижимому случаю, позволяющему характеризовать минимальную амплитуду простым соотношением. На практике минимальная амплитуда будет находиться в диапазоне , поэтому использование PL является вполне обоснованным.

Заключение

Предложенный алгоритм обнаружения гармонических составляющих позволяет определить частоту гармонической помехи с целью её дальнейшей фильтрации. Как теоретический анализ, так и имитационное моделирование, показали, что с помощью предлагаемого алгоритма, частота сигнала в большинстве случаев может быть определена с абсолютной погрешностью f<0.5fs. Приведенный расчет предельной погрешности может служить основой для выбора спектрального окна, наилучшим образом подходящего для обнаружения гармонических составляющих в интересующем исследователя диапазоне амплитуд. В качестве универсального окна, позволяющего добиться приемлемого результата в широком диапазоне амплитуд (от 0 до 40 дБ) можно рекомендовать окно Ханна.

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

Литература

1. Айфичер Э.С., Джервис Б.У. Цифровая обработка сигналов: практический подход, 2-е издание: Пер. с англ. - М.: Издательский дом «Вильямс», 2004. 992 c.

2. Бендат Дж., Пирсел А. Прикладной анализ случайных данных: Пер. с англ. - М.: Мир, 1989. 540 с.

3. Harris F.J. On the use of windows for harmonic analysis with the discrete Fourier transform, Proc. IEEE, 66(1), 51-84, 1978.

4. Oppenheim A.V., and R.W. Schafer, Discrete-Time Signal Processing, Prentice-Hall, 1999.

Размещено на Allbest.ru