ПРИМЕНЕНИЕ РАЗРЕЖЕННОГО ПРЕОБРАЗОВАНИЯ ФУРЬЕ ДЛЯ АНАЛИЗА УЗКОПОЛОСНЫХ РАДИОСИГНАЛОВ Текст научной статьи по специальности «Электротехника, электронная техника, информационные технологии»

www.mai.ru/science/trudy/

Труды МАИ. Выпуск №84

УДК: 621.396.1

Применение разреженного преобразования Фурье для анализа

узкополосных радиосигналов

Трофимов Д.В.*, Баев А.Б.**

Московский авиационный институт (национальный исследовательский университет), МАИ, Волоколамское шоссе, 4, Москва, A-80, ГСП-3, 125993, Россия

*e-mail: radiovtuz@bk.ru **e-mail: baev@mai.ru

Аннотация

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

Ключевые слова: разреженность, спектральный анализ, преобразование Фурье, прореживание, мониторинг электромагнитных излучений.

Введение

В процессе развития цифровой техники из-за увеличения объёма обрабатываемых данных возникла необходимость в разработке таких алгоритмов преобразования Фурье (ПФ), которые требуют меньших вычислительных затрат по сравнению с обычным алгоритмом быстрого преобразования Фурье (БПФ). Если спектр сигнала разрежен, то есть, содержит небольшое число значимых гармонических составляющих, то такой сигнал можно аппроксимировать, вычисляя только наиболее значимые коэффициенты ряда Фурье [1]. Для этого используются алгоритмы разреженного преобразования Фурье. Некоторые, наиболее простые алгоритмы РПФ, позволяющие определять частоту комплексного гармонического сигнала, кратко рассмотрены в [1]. Определение частоты одной гармоники - это узкая задача, в большинстве же практических задач требуется быстро найти частоты, амплитуды и фазы большого числа гармоник. Одной из таких задач является контроль электромагнитных излучений с целью обеспечения электромагнитной совместимости бортовой радиоэлектронной аппаратуры.

Алгоритм с прореживанием во временной области

Подробное описание алгоритма РПФ с примером его аппаратной реализации представлено в [2]. Выигрыш в объёме вычислений достигается за счёт того, что число отсчетов полученного спектра существенно меньше, чем длительность анализируемого сигнала, так как БПФ берётся от прореженного во временной области сигнала. Особенностью алгоритма является то, что требуется заранее задать предполагаемое число значимых гармоник К, а также то, что исследуемый сигнал

должен быть комплексным. При прореживании действительного сигнала, обладающего двухсторонним спектром, могут возникнуть ошибки из-за наложения спектров. Например, если прореживать синусоидальный сигнал с шагом, кратным половине периода, в частотной области будет производиться сложение отсчётов нечётно-симметричного спектра (мнимой части спектра). Эти отсчеты компенсируют друг друга, и амплитудный спектр получается равным нулю -гармоника исчезает. В случае с прореживанием косинуса можно по той же причине получить удвоенную амплитуду.

Введём следующие обозначения: s[n] - исследуемый сигнал; S[m] - спектр исследуемого сигнала; S"[m] - аппроксимация спектра исследуемого сигнала, полученная в результате РПФ; N - длительность сигнала и спектра, n = [0...N-1], m = [0...N-1]. Длительность сигнала N должна удовлетворять условию: она должна быть такой, чтобы её можно было представить в виде произведения Q взаимно простых чисел N1 ...Nq.

Q

N = П N. (1)

q=l

Зададимся сигналом, являющимся копией сигнала s[n], задержанной на t отсчётов:

st [n] = s[n +1 ]. (2)

Далее исследуемый сигнал и его задержанная копия (2) прореживаются в Lq раз:

Lq = N / Nq , (3)

где q = [1... Q]. Процесс прореживания можно записать так:

ия [пч ] = 8[п ■ 1Ч ], V[пя] = [п ■ ЬЧ ].

Всего после прореживания получается 2•Q новых сигналов: uq[nq] и vq[vq], где пС1 = [0...^-1], q = [1..^], Nq - длительность q-го сигнала. Надо заметить, что Q сигналов происходят из незадержанного сигнала s[n], и Q сигналов из задержанного сигнала ^[я]. От каждого прореженного сигнала вычисляется дискретное преобразование Фурье (ДПФ):

1 N-1 Тп"т"

и [т ] = — У и [п ] ■ е

ЯL ЯJ ЛТ I ЯL Я-1

N пЯ=о

Я ' 2, (5)

V т ] = — У V п ] ■ е ^

^ Ял\т ЯЯ -

N п„=0

Я

"Я

где шС1 = [0...^-1]. В иностранной литературе отсчёты спектров (5) для удобства принято называть «корзинами». В одном из спектров (5), например в Ц^ш^, определяются номера гармоник с амплитудой большей, чем пороговый уровень А0, заданный заранее. Из этих номеров составляется матрица I. Эта матрица содержит Q строк, а количество столбцов Н определяется количеством значимых гармоник. Максимально возможное количество столбцов будет равно наибольшей длительности спектра ^[ш^, то есть равно наибольшему из взаимно простых чисел: Нтах = тах{^}. Если в q-м спектре все отсчёты меньше А0, то вся q-тая строка матрицы остаётся «пустой». Далее происходит последовательное сравнение амплитуд отсчётов спектров (5):

и [ ]| - К [ ]| (6)

где к = [1...Я], е - некоторая заранее заданная величина. В случае идеального сигнала без шумов следует принять е = 0. Если условие (6) выполняется, то в новую матрицу J добавляется элемент, равный 1ф

н = н. (7)

Количество «непустых» элементов в матрице J, как правило, меньше чем, в I (так как при наличии шумов условие (6) не всегда выполняется), но её максимальный размер будет такой же, как у I. Каждому элементу Jqк матицы J соответствует некоторая частота в спектре 5[т]. Они связаны формулой:

[ .1Чн н ]/ V [ 3ЧН Н ])

2ж-Г

ХЧ,Н = N 4 4 , (8)

где - элемент матрицы X, в которую записываются все найденные частоты. Одна и та же частота может быть обнаружена несколько раз, то есть в матрице X могут быть одинаковые элементы. Для того чтобы составить матрицу М (она содержит три строки и неизвестное заранее число столбцов), где хранятся окончательные результаты работы алгоритма, нужно записать в неё только разные частоты из матрицы X. То есть в первой строке матрицы М не должно быть повторяющихся элементов. Остальные две строки используются для записи в них комплексных амплитуд гармоник с обнаруженными частотами (для задержанного и не задержанного сигнала). Обобщая сказанное можно записать так:

Ми = Хчн; МХг = ид];М3, = Уд], (9)

где I - текущий номер обнаруженной гармоники, при условии, что не повторяются. На этом первый цикл вычислений заканчивается. Далее оценивается количество столбцов в матрице М. Если оно больше или равно К, то вычисления прекращаются. А если нет, то происходит вычитание из сигнала и его задержанной копии гармоник с найденными частотами, амплитудами и фазами. Затем цикл вычислений повторяется.

Особенностью данного алгоритма является то, что он способен накапливать ошибки. Это происходит следующим образом: если частота гармоники была неправильно определена, то при вычитании её из сигнала вместо компенсации происходит добавление в сигнал новой гармоники. Ошибки возникают при наличии в сигнале шумов. Как показало моделирование, алгоритм устойчиво работает только при отношениях сигнал/шум более 50 дБ. Алгоритм можно улучшить, если ввести в него дополнительную проверку найденных частот. Выражение (8) позволяет определить истинную частоту по фазовому сдвигу (по отношению корзин спектров (5)), но сам номер корзины так же указывает на возможные значения частоты. Из-за наложения спектра корзине с номером Мь могут соответствовать несколько частот в спектре исходного сигнала:

М = М., М + N , М + 2 N ,...М + N ■

воз Ь? Ь я Ь я Ь Я

с / Л \

N -1

целое

N

V V я У

(10)

Воспользоваться этим свойством можно следующим образом: по формуле (1 0) составляется вектор Мвощ[г] длиной ^ = 1+^-(целое(М^) - 1), г = [1...Лу. Далее находится разница между каждым его элементом и значением частоты, найденным по (8):

MA[r] = Meo3[r] - X, h . (11)

Определяется минимальный элемент (12) и соответствующий ему номер rmin:

Ma min[rmin] = mm{Ma }. (12)

За истинную оценку частоты Xh принимается элемент вектора Мвоз с номером rmin:

X „ = M [r ]. (13)

q, h воз L min J V /

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

Алгоритм с прореживанием в частотной области Другим подходом к уменьшению количества вычислений при нахождении спектра сигналов является прореживание спектра [3]. Данный подход основан на дробном преобразовании Фурье и изначально применялся для обработки импульсов с линейной частотной модуляцией (ЛЧМ-импульсов) в радиолокации [4]. Исследуемый сигнал может быть действительным или комплексным. В соответствие со свойствами преобразования Фурье, прореживанию в частотной области соответствует наложение сигнала во временной области. Введём следующие обозначения: z[b] - сигнал, получаемый из исходного, путём наложения; Z[v] - спектр сигнала z[b]; b = [0...B - 1]; v = [0...B - 1], где B - длительность сигнала z[b]. Спектр S[m] связан со спектром Z[v] следующим образом:

N'

Z[v] = S

v — B

= S[v • L], (14)

где Ь=Ы1В - коэффициент прореживания спектра (Ы должно быть кратно В). Во временной области спектру, определённому в (14), соответствует сигнал 2[Ь],

полученный как сумма Ь периодов сигнала s[n], задержанных на В отсчётов относительно друг друга:

1 ь-1

Ф] = - ■У s[Ь +1 ■ Ь]. (15)

Ь 1=0

Спектр 7\у] можно найти с помощью ДПФ сигнала (15):

1 в-1 -1 ,Ь

= - У 2[Ь] ■ е (16)

В Ь=о

Таким образом, спектр вычисляется для сигнала длиной В, а не N, что уменьшает объём вычислений. После нахождения прореженного спектра по формуле (16), нужно сравнить амплитуды найденных гармоник с пороговым уровнем А0. Номера отсчётов с амплитудой больше А0 сохраняются в векторе 0[р], р = [1...Р], где Р - количество значимых гармоник. Далее нужно восстановить истинные номера М1...МР гармоник в спектре £[ш], т = [0...У-1]. Для этого используется функция соответствия И[п], заданная как:

Н[п] =

в ■ п п

_ N _ _ Ь _

(17)

Функция к[п] = [0...В -1], п = [0...У -1], позволяет установить связь между отсчётами спектра 7\у] и 5[ш]. Она имеет ступенчатый вид, причём высота ступенек растёт дискретно с шагом единица, а ширина каждой ступеньки равна Ь. Главное свойство функции И[п] в том, что нескольким значениям аргумента соответствует одно значение функции. Истинный номер гармоники определяется следующим образом: если значение Мр элемента с номером р из вектора G[p] совпадает с одним из значений функции И[п], то всем элементам вектора 5"[ш], попадающим под Мр-ю

«ступеньку» присваивается значение амплитуды отсчёта Z[Mk]. Каждой гармонике из спектра 8[т] соответствует Ь гармоник одинаковой амплитуды в спектре 8'[т]:

S[M ].. Б[М + L] = Z[Н[M ] ].. Z[Н[M + L]] (18)

Только крайний отсчёт с наименьшим номером будет соответствовать истинному отсчёту в спектре 8[т]. Рассмотренный алгоритм обладает очевидным недостатком: он не способен обнаруживать гармоники с номерами не кратными Ь. Для устранения этой проблемы можно искусственно расширить спектр. Проще всего свернуть в частотной области каждый отсчёт (дельта-функцию) с прямоугольным спектром Щт]

Ж [т] = гваь [т] (19)

оконной функции w[n] вида:

МШ] = Ь • sinc(лnL). (20)

Во временной области этому соответствует умножение сигнала на оконную функцию, в результате чего получается новый сигнал с[п]:

с[п] = ь[п] • м\п] ^ С[т] = Б[т] * Ж[п]. (21)

При прореживании спектра С[т], даже если номер гармоники М не равен у-Ь, в любом случае останется ненулевой отсчёт. При этом должно соблюдаться условие, что промежуток между двумя соседними отсчётами в спектре исходного сигнала не должен быть меньше, чем Ь. Иначе два прямоугольника «наложатся» друг на друга, что приведёт к ошибкам. В определённых случаях можно устранить и это препятствие. Для этого используется определённая перестановка элементов вектора s[n], определяемая формулой:

\n] = s\(n • a) mod N] = s\na ], (22)

где sa[n] - новый вектор, состоящий из элементов вектора s[n], а - параметр перестановки, a = \1..N ], n = \0..N -1], na = \1...N" ]. Число na является сравнимым по модулю с числом (n•a), то есть, при делении на N они должны давать одинаковый остаток. Понятно, что при одном и том же (n • a) таких чисел na будет бесконечно много, поэтому надо выбрать минимальное из них. А минимальным числом na будет остаток от деления (n • a) на N:

'n a

nCT = остаток

(23)

v N

N

Согласно (23) производится прореживание сигнала s[n] в о раз, в результате чего получается сигнал ха[и], а затем такие сигналы составляются друг за другом о раз. В итоге частота всех гармоник в спектре увеличивается в о раз, спектр как бы «растягивается». При этом желательно выбирать о таким, чтобы N делилось на него без остатка, иначе может возникнуть растекание спектра.

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

ha\n] =

a • B • n a • n

_ N _ _ L _

(24)

Если сравнить (17) и (24), то можно выявить следующие отличия: • функция ^ [п] растёт быстрее;

К[п] может принимать значение больше, чем В. Отсюда следует, что частота гармоник в анализируемом сигнале не должна превышать значения N/0;

"a

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

Применение алгоритмов РПФ для исследования модели сигнала

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

Таблица 1.

Параметр Для алгоритма с прореживанием по времени Для алгоритма с прореживанием по частоте

Тип сигнала Комплексный

Длительность сигнала N 101102 = 10302 104

Частоты гармоник 100, 200, 210, 380

Амплитуды гармоник 2, 1, 3, 5, 4

Фазы гармоник п/3, п/6, п/2, п, п/5

Шум Белый гауссовский, отношение сигнал/шум 20 дБ

Результат вычисления спектра для алгоритма с прореживанием по времени представлен на рис. 1. Представленные результаты демонстрируют влияние

задержки сигнала (параметр ?) и применение проверки найденных частот (13). На рис. 1а) и 1в) задержка составляет один отсчёт, на рис. 1б) и 1г) - десять отсчётов. На рис. 1 а и 1 б) спектр определялся без поиска ближайшей из возможных частот, а на рис. 1 в) и 1 г) - с поиском.

а)

в)

б)

—© РПФ —»< ДПФ 1 к

..... \ Т

ч »

1

< > А

_ .....

_]_I_I_I_I_1_

Г)

О 50 100 150 200 256 300 350 400 450 500

Номер отсчёта, т

Рис. 1. Спектры модели сигнала, найденные при помощи алгоритма РПФ с

прореживанием по времени.

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

так же существенно улучшает результат. Наилучший результат получается при сочетании обоих способов.

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

а)

в)

г)

Рис. 2. Спектры модели сигнала, найденные при помощи алгоритма РПФ

с прореживанием по частоте.

На рис. 2а) и 2в) количество интервалов В = 100, на рис. 2б) и 2г) В = 1000. На рис. 2а) и 2б) не используется перестановка элементов исследуемого сигнала. На

рис. 2в) и 2г) параметр перестановки равен 10, то есть спектр исходного сигнала «раздвигается» в 10 раз. Спектр, определённый на рис. 2а, практически не пригоден для анализа сигнала. Указаны лишь широкие области, где могут быть гармоники. Гармоники с частотами 200 и 210 вообще не различимы. Точность определения диапазона частот, где находятся гармоники, может быть улучшена как за счет увеличения количества интервалов до 1000 (рис. 2в), так и путем введения перестановки элементов сигнала (рис. 2б). Но второй способ более предпочтителен, так как не увеличивается размер вычисляемого ДПФ. Применить такую перестановку можно, так как номера гармоник в исследуемом сигнале существенно меньше длительности сигнала. При увеличении частот всех гармоник наибольшая частота не должна иметь номер больше N. В данном случае наибольшее значение параметра а равно ММтах = 104/3 80 « 26. Тем не менее, алгоритм указывает лишь интервалы частот (конечно, более узкие по сравнению с рис. 2а), где могли быть исходные гармоники. Как и в случае с алгоритмом с прореживанием по времени, наилучшее разрешение по частоте обеспечивается при сочетании двух способов. Аппроксимация спектра на рис. 2г) ближе всего к истинному спектру сигнала.

Применение алгоритмов разреженного преобразования Фурье для анализа

реальных сигналов

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

данных в персональном компьютере. На рис. 3 представлены аппроксимации спектра, полученные при помощи алгоритма с прореживанием по времени (рис. 3 а) и алгоритма с прореживанием по частоте (рис. 3б). На рис. 4 показаны увеличенные фрагменты спектра вблизи частоты 60 МГц для двух алгоритмов. Частота гармоники, равная 60 МГц, с помощью первого алгоритма РПФ определена как

61 МГц, то есть погрешность составляет примерно 2 %. При этом точно определена её амплитуда - 87 мВ. Так же было обнаружено 12 несуществующих гармоник на разных частотах и с разными амплитудами. Алгоритм с прореживанием по частоте (рис. 3б) не даёт непредсказуемых ошибок, но принципиально не способен точно указать на истинную частоту. Из рис. 4б) следует только, что в диапазоне 60 -

62 МГц находится гармоника с амплитудой 73 мВ. Погрешность определения амплитуды при этом составляет 16 %.

а)

Рис. 3. Оценка спектра измеренного сигнала, полученная с помощью алгоритма с

прореживанием по времени а) и алгоритма с прореживанием по частоте б).

Рис. 4. Фрагменты спектра измеренного сигнала, полученные с помощью алгоритма с прореживанием по времени а) и алгоритма с прореживанием по частоте б).

Выводы

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

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

Выбор алгоритма РПФ так же зависит от решаемой задачи. Алгоритм с прореживанием по времени и задержкой способен находить много гармоник, причём разрешение по частоте принципиально не ограничено. При этом он обладает низкой устойчивостью к шумам и чувствителен к растеканию спектра. Он может быть использован там, где уровень шума очень низкий. Неудобство состоит и в том, что исследуемый сигнал должен быть комплексным или состоять из двух квадратурных составляющих, которые могут быть объединены в комплексный аналитический сигнал без существенного увеличения вычислительных затрат. Применение такого алгоритма к действительным сигналам неэффективно, поскольку требует определения преобразования Гильберта, которое в свою очередь приводит к увеличению объёма вычислений. Низкая устойчивость к шумам алгоритма объясняется тем, что он способен накапливать ошибки. Это происходит из-за того, что невозможно скомпенсировать гармонику, вычитая из неё гармонику с другой частотой. Введение в алгоритм поправки за счёт определения ближайшей из возможных частот повышает его устойчивость к шумам, но несколько увеличивает время работы. Результат работы будет зависеть и от того, на какие взаимно простые числа будет разделена длительность сигнала Я, а также от заданного заранее предполагаемого числа значимых гармоник К. Если К выбрано слишком маленьким, он может не найти все гармоники, а если слишком большое - время работы не оправдано возрастает.

Алгоритм с прореживанием по частоте использует все отсчёты сигнала. Но, несмотря на это, выигрыш в скорости обработки достигается за счёт уменьшения размера ДПФ или БПФ. При этом ухудшается разрешение по частоте во столько раз, во сколько был уменьшен размер БПФ (в L раз). Выбор разрешения по частоте зависит от разреженности сигнала и требуемой точности. При повышении точности снижается скорость работы. Разрешение по частоте можно улучшить, вводя перестановку элементов сигнала, но при этом нужно следить, чтобы не нарушалась теорема Котельникова, так как частота всех гармоник повышается в о раз. Алгоритм способен обрабатывать как комплексные, так и действительные сигналы, что удобно при применении к реальным сигналам. Данный метод может быть эффективно применен, например, при сканировании эфира, когда возникает необходимость быстро узнать, в каком диапазоне частот появилось излучение.

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

Библиографический список

1. Gilbert A.C., et al., Recent Developments in the Sparse Fourier Transform, Signal Processing Magazine, 2014, vol. 31, №. 5, pp. 91-100.

2. H. Hassanieh, et al., "GHz-Wide Sensing and Decoding Using the Sparse Fourier Transform", IEEE International Conference on Computer Communications (INFOCOM), Toronto Canada, April 2014, pp. 2256-2264.

3. Shengheng Liu, et al., "Sparse Discrete Fractional Fourier Transform and Its Applications", IEEE Transactions on Signal Processing, 2014, vol. 62, № 24, pp. 65826595.

4. Soo-Chang Pei, et al., "A New Discrete Fractional Fourier Transform

Based on Constrained Eigendecomposition of DFT Matrix by Largrange Multiplier Method", IEEE Transactions on Circuits and Systems - II: Analog and Digital Signal Processing, September 1999, vol. 46, № 9, pp. 1240-1245.