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

Радиолокация и радионавигация

УДК 621.396.96

А. В. Бархатов, А. С. Козлов Санкт-Петербургский государственный электротехнический университет "ЛЭТИ" им. В. И. Ульянова (Ленина)

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

Описан алгоритм быстрого вычисления частотно-временной функции в радиолокационной станции (РЛС). Сделан акцент на особенностях программной реализации алгоритма. Описаны особенности параллельных вычислений функции на процессорах игровых видеокарт. Определено время, затрачиваемое на параллельные вычисления 16 частотно-временных функций в многоканальной РЛС.

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

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

N-1

d)) =

X * (n) 4f (n -1) e-j2wdn'N

n = 0

, (1)

где l - дискретизированная задержка (lag); d -дискретизированный частотный сдвиг; N - количество отсчетов принятого сигнала в обрабатыва-

емом блоке; n = 0, N -1 - номер отсчета сигнала; s (n) - принятый сигнал; sref (n -1) - цифровой

опорный (referenced) сигнал; - символ комплексного сопряжения.

Дискретизированная задержка l представляет собой задержку сигнала т, дискретизированную на интервал дискретизации (sampling) Ts : l = round (т/Ts)

(round (•) - операция взятия целой части аргумента). Дискретизированный частотный сдвиг аналогично связан с доплеровским смещением /D и длительностью анализируемой выборки сигнала NTs :

d = round (fD NTs) = round [(/D/ fs )N ], где / = 1/Ts - частота дискретизации.

Если s(n) совпадает с опорным сигналом sref (n -1), то (1) является оценкой функции неопределенности (ФН) сигнала sref (n -1) на интервале задержек 0...LTs и частот - D/s/ N... D/s/ N. На практике в принятом сигнале обычно содержится множество копий sref (n -1), сдвинутых

по задержке и по частоте.

В зарубежной литературе функцию (1) называют "поверхностью амплитуда-дальность-Доплер" (amplitude-range-Doppler surface) или просто "поверхностью дальность-Доплер" из-за прямой связи задержки с дальностью и обусловленностью частотного сдвига эффектом Доплера. В отечественных источниках ее называют корреляционно-фильтровой функцией, взаимной функцией неопределенности, двумерной взаимно корреляционной функцией, функцией отклика и т. д. Далее будем называть (1) частотно-временной функцией (ЧВФ). На рис. 1 приведен пример ЧВФ, нормированной (normalized) на свое максимальное значение, выраженной в децибелах.

Алгоритм расчета частотно-временной функции. Различные способы расчета ЧВФ, учитыва-

В публикации использованы результаты ОКР "Разработка пассивного когерентного локационного комплекса для охраны важных объектов", выполняемой СПбГЭТУ "ЛЭТИ" по договору с АО «НИИ "Вектор"» в рамках комплексного проекта по созданию высокотехнологичного производства при финансовой поддержке работ по проекту Министерством образования и науки Российской Федерации (постановление Правительства Российской Федерации от 9 апреля 2010 г. № 218).

42

© Бархатов А. В., Козлов А. С., 2015

I* п

Щй 0\ -

/в, Гц 0 Рис. 1

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

Пусть принятый сигнал 5 содержит N отсчетов, а опорный сигнал 5ге£ - N + Ь отсчетов. На

рис. 2, а блоки отсчетов условно показаны прямоугольниками.

Блок из N отсчетов сигнала 5 разбивается на М сегментов по N отсчетов в каждом. Аналогично разбиваются N отсчетов сигнала 5ге£ и комплексной экспоненты. При этом можно записать:

М-1N'-1 d)|= £ £ 5 ' + П')х т=0 п '=0

х 5^ (тМ' + п' -1)е- ]2Ы(т^ +П)(Ш'|, (2)

где т = 0, ..., М -1 - номер сегмента; п - номер отсчета внутри сегмента.

N отсчетов

5 N отсчетов N' отсчетов | N' отсчетов

Сегмент 0 N + Ь Сегмент | М -1 отсчетов |

г-— | —■>

Ь отсчетов N отсчетов N' отсчетов ! W отсчетов

| т-й сегмент >

5ге£

Нули N отсчетов Нули

т-й сегмент

Ь отсчетов N отсчетов Нули

б

Рис. 2

С учетом того, что

- ] 2%d т'+п ' )/(ш') = = -]2%dmn'/(Мn) -]2я<аП'¡(ш') =

= - ] 2^т/м е-] '/(mn' ) — е е ,

получим

, d ) =

т=0 _

х ,Гет Ш' + п' -1)е-]2ЫП'!(Ш'}]

М -1

£

т=0

N'-1

е-]'2™1т1М £ 5 ^ ' + П)х п '=0

Рассмотрим случай Б ^ М, что означает, что максимальное значение модуля частотного сдвига много меньше частоты повторения сегментов /м :

, Б } ( , / М Л •^тах = J < I/М = N ' ).

При этом максимальное значение модуля показателя экспоненты в (2)

[2%сП'/(МЫ')]тах = 2%0 (N' - 1)/(Ш')« 2я£/М « 0.

Т - ][2Ып ' /(MN' )]тях - ]0 л П

Тогда е ■1тах « е 3 = 1. При мень-

ших значениях d и п' это равенство тем более выполняется.

В этом случае ЧВФ

И, d)| =

М-1

£

т=0

N -1

е~]2%Ап/М £ 5 ^' + п )х П '=0

х 5*е£ (тИ' + п -1)]

(3)

Обозначим

N'-1

гт (I) = £ 5 (mN' + п') 5*е£ (mN' + П -1) (4)

п =0

- дискретную линейную свертку т-го сегмента 5 и т-го сегмента 5*^ (рис. 2, б), вычисленную в ограниченном диапазоне сдвигов 0 < I < Ь.

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

а

5

последовательность образована отчетами в т-м

сегменте 5, вторая - отсчетами в т-м сегменте

*

5ге£. Сокращение вычислительных затрат достигается за счет вычисления ДПФ с помощью алгоритма БПФ. Другими словами, для быстрого вычисления (4) нужно хотя бы один из сегментов (или принятого, или опорного сигналов) дополнить нулевыми отсчетами, выполнить БПФ для каждого из сегментов, перемножить полученные спектры и выполнить обратное БПФ. При программной реализации массивы с отсчетами сигналов должны иметь одинаковую длину, поэтому нулевыми отсчетами дополняются оба сегмента. При этом, во-первых, сегмент принятого сигнала 5 сначала дополняется Ь нулевыми отсчетами слева (см. рис. 2, б), во-вторых, количество нулевых отсчетов, которыми дополняются оба сегмента справа, выбирается таким, чтобы длина дополненных массивов не только равнялась или превышала Ь + N' + N' — 1 (чтобы круговая свертка совпадала с линейной), но и была степенью двойки - при этом обеспечивается наилучшая вычислительная эффективность БПФ.

С учетом (4) перепишем (3):

И, <01 = М——Ъ (I) е~^т1М . (5)

т=0

Для каждого значения задержки I (5) представляет собой ДПФ вектора

[г0 (I), 1 (I), •••, М—1 (I)]т

( т - символ транспонирования), которое также можно рассчитать по алгоритму БПФ для сокращения вычислений. Таким образом, для каждого значения I сечение ЧВФ вычисляется с помощью БПФ. На рис. 3 проиллюстрировано, как совпадающие по задержке отсчеты гт (1) объединяются в вектор для последующего расчета 1-го сечения ЧВФ.

Перед БПФ вектор [г0(1), 1\(1), •.., гМ—1 (1)]т должен быть дополнен нулевыми отсчетами так,

r0 (0) r0 (1) r0(L)

r1 (0) r (1) r (L)

rm (0) rm (D rm (L )

Гм-1 (0) rM-1 (1) rM-1 (L )

Рис. 3

чтобы длина дополненного вектора, во-первых, была не меньше 2M -1 и, во-вторых, была степенью двойки. Для снижения уровня боковых лепестков ЧВФ перед дополнением нулевыми отсчетами вектор следует умножить на оконную функцию (например, Хемминга) длиной M.

На основании изложенного алгоритм вычисления ЧВФ включает следующие операции:

1. Получение векторов s длиной N и sref длиной L + N.

2. Комплексное сопряжение вектора sref.

3. Разделение вектора s длиной N на M векторов sm длиной N' (см. рис. 2).

*

4. Разделение вектора sref на M векторов

* I

sref m длиной N + L (см. рис. 2).

5. Определение числа точек БПФ rm (l) (4):

NFFTr = 2nextpow2(2N'+L-1) (FFT - Fast Fourier Transformation - быстрое преобразование Фурье; nextpow2(n) - функция, возвращающая минимальное p, такое, что 2p > n).

6. Дополнение каждого вектора sm L нулями слева и NFFTr - L - N' нулями справа.

*

7. Дополнение каждого вектора sref m NFFTr - L - N' нулями справа.

*

8. БПФ для всех векторов sm и sref m (далее -БПФ-8).

9. Перемножение векторов-результатов БПФ-8 для совпадающих m.

10. Обратное БПФ (далее - БПФ-10) для всех M векторов-результатов перемножения из п. 9.

11. Отбрасывание в каждом векторе-результате БПФ-10 последних NFFTr -(L +1) отсчетов.

12. Формирование матрицы Rmi с размерами M х (L +1) из векторов-результатов БПФ-10.

13. Умножение каждого столбца матрицы Rmi на оконную функцию длиной M.

14. Определение количества точек БПФ

Y(l, d) (5): Nffty = 2nextpow2(2M

15. Дополнение всех столбцов матрицы Rmi Nfft y - M нулями.

16. БПФ (далее - БПФ-16) для каждого столбца матрицы Rmi с размерами Nfft y х(L +1).

17. Формирование новой матрицы с размерами (2D + 1)х( L +1).

18. Копирование в первые D столбцов сформированной матрицы из п. 17 последних D столбцов матрицы-результата из п. 16, а в последние D +1 столбцов указанной матрицы - первых D +1 столбцов матрицы-результата из п. 16.

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

Видеокарты привлекательны тем, что их архитектура позволяет организовать параллельные вычисления за счет многопроцессорности. Например, используемая в описанной далее реализации относительно недорогая видеокарта ASUS Radeon R9 290 содержит 2560 универсальных вычислительных процессоров. Кроме того, при параллельных вычислениях для увеличения суммарного объема памяти или для повышения производительности можно использовать не одну, а две и более видеокарт.

Первым этапом вычисления ЧВФ на видеокарте является подготовка данных для параллельных вычислений БПФ-8. Подготовка заключается в копировании отсчетов принятого и опорного сигналов в память видеокарты, разделении их на сегменты и дополнении нулями. Для минимизации количества операций по работе с памятью видеокарты (поскольку они являются самыми медленными) следует не копировать отсчеты в память и затем перемещать их там, дополняя нулями (при этом на каждый элемент обрабатываемого массива придется как минимум две операции копирования), а предварительно заполнить память нулями и, используя простой алгоритм, скопировать каждый элемент всего один раз. На рис. 4 это упрощенно проиллюстрировано для сигнала 5 (здесь подразумевается, что данные в память видеокарты копируются из оперативной памяти компьютера).

Несмотря на то, что при таком варианте приходится копировать каждый сегмент по отдельности (в отличие от копирования целиком всего блока от-

Сегменты зт

в памяти видеокарты перед БПФ-8

Сигнал ^

в оперативной__

памяти * у ; * ' N N ячеек ячеек Рис. 4

счетов сигнала), временные потери оказываются невелики, так как передача данных осуществляется по технологии прямого доступа к памяти.

В рассмотренной реализации компьютерная программа вычисления ЧВФ создана с помощью ОрепСЬ (фреймворк для написания компьютер -ных программ, связанных с параллельными вычислениями, содержащий свои язык программирования и интерфейс программирования приложений). БПФ (шаги 8 и 16 алгоритма вычисления ЧВФ), а также обратное БПФ (шаг 10 алгоритма) осуществлялись с помощью библиотеки сШБТ ОрепСЬ (при этом настройка БПФ сводилась к указанию количества точек, типа данных, буферов для хранения входных и выходных данных). Требовалось обеспечить параллельное вычисление ЧВФ для 16 приемных каналов. Количество необходимой для обработки 16 каналов памяти превышало объем памяти одной видеокарты, поэтому параллельные вычисления выполнялись на двух видеокартах. При этом пришлось создать 2 контекста (рабочих пространства) для библиотеки сШБТ, каждый из которых содержал по одной видеокарте. На скорости вычислений использование двух контекстов вместо одного не сказывалось.

Перемножение спектров (шаг 9 алгоритма) выполнялось ядром ОрепСЬ. ОрепСЬ не содержит операций для работы с комплексными числами, поэтому они были реализованы с помощью специально написанных тИпе-функций. Для представления комплексного числа в ОрепСЬ использован тип данных Аоа12, содержащий два числа с плавающей точкой.

После выполнения обратного БПФ-10 выходной буфер данных (здесь и далее - для каждого приемного канала) содержал М векторов по ^рртг комплексных отсчетов в каждом. В процессе реализации шагов 11-15 алгоритма эти векторы преобразовывались в матрицу Ят1 с размерами Т ^ +1).

Отсчеты матрицы содержались во входном буфере перед выполнением БПФ-16. В рассмотренной реализации оба буфера (выходной буфер после

wFFT ячеек "fft ячеек

V ч

0 0 0 0 0 0 0 0

/

обратного БПФ-10 и входной буфер перед БПФ-16) не могли быть размещены в памяти видеокарты одновременно. В связи с этим создавался промежуточный буфер меньшего размера, в который копировались Ь +1 отсчетов каждого из М векторов, после чего выходной буфер обратного БПФ-10 уничтожался. Затем создавался и заполнялся нулями буфер для матрицы В.т1, после чего отсчеты из промежуточного буфера копировались в новый буфер по аналогии с тем, как это делалось при загрузке данных в память видеокарты перед выполнением БПФ-8 (см. рис. 4).

Матрицу Ят[ следует ориентировать в памяти видеокарты так, как показано на иллюстративном примере (рис. 5) для матрицы размером 3 х 4, чтобы друг за другом располагались отсчеты одной строки (рис. 5, а), а не столбца (рис. 5, б), несмотря на то, что БПФ выполняется для столбцов матрицы. Заливкой на рис. 5 выделены ячейки с элементами первой строки. Для операции выполнения БПФ расположение элементов в памяти не играет роли, но при последующем извлечении данных (в том числе из памяти видеокарты) для выполнения шагов 17, 18 алгоритма это важно, и размещение по строкам заметно экономит время.

r11 r12 r13 r14 r21 r22 r23 r24 r31 r32 r33 r34

а

r11 r21 r31 r12 r22 r32 r13 r23 r33 r14 r24 r34

б

Рис. 5

Время вычисления частотно-временной функции. Описанным способом выполнялось параллельное вычисление 16 ЧВФ. Блок комплексных отсчетов для каждого из 16 принятых сигналов 5 имел длину N = 2 308 198 отсчетов. Опорный

сигнал 5ге£ содержал N + Ь = 2 308 198 + 4 572 = = 2 312 770 отсчетов. Каждый сигнал разделялся на М = 1262 сегментов по N' = 1829 комплексных отсчетов в каждом (перекрывающиеся сегменты опорного сигнала содержали по N + Ь = 1829 + 4572 = 6401 отсчету). Количество

точек БПФ-8 = 16 384. На шаге 8 алгорит-

ма для всех сигналов (16 принятых и одного опорного) суммарно выполнялось 21 454 БПФ на 16 384 точках. На шаге 10 алгоритма для всех сигналов выполнялось 20 192 обратных БПФ на 16 384 точках. Количество точек БПФ-16 =

= 4 096 (матрица Rmi на этом шаге алгоритма имела размеры 4096 х 4573). На шаге 16 алгоритма для всех приемных каналов суммарно выполнялось 73 168 БПФ на 4 096 точках. Размеры выходной матрицы, представляющей ЧВФ, 1537 х 4573. Расчет 16 ЧВФ на двух видеокартах ASUS Radeon R9 290 с учетом загрузки данных в память видеокарт и извлечения результатов в оперативную память компьютера занял 660 мс.

Полученный результат - время, затрачиваемое на расчет ЧВФ, позволяет, на взгляд авторов, рассматривать описаный вариант аппаратно-программной реализации цифровой обработки как приемлемый для использования в реальной практике. В частности, при интервале дискретизации порядка 100 нс блок отсчетов сигнала длиной, аналогичной приведенной ранее (несколько миллионов отсчетов), соответствует времени накопления несколько сот миллисекунд, соизмеримому со временем обработки на видеокартах. Можно сформулировать по-другому. В рассмотренном примере две видеокарты обработали около 40 млн отсчетов за 660 мс. Поскольку каждый отсчет является 64-разрядным, обрабатываемый поток составил 3.8 Гбит/с, что также является неплохим показателем. Рассмотренный вариант реализации вычисления ЧВФ позволяет увеличивать количество обрабатываемых приемных каналов без повышения времени вычислений за счет добавления видеокарт. С точки зрения стоимости рассмотренный вариант также привлекателен, так как обычно специализированный вычислитель стоит существенно дороже видеокарт. Программный код, написанный в OpenCL, не привязан к конкретной аппаратуре, что позволяет использовать видеокарты различных производителей. Во многих случаях гражданского применения, когда использование на РЛС компьютеров допускается техническими требованиями, применение видеокарт для цифровой обработки сигналов является реальной альтернативой.

СПИСОК ЛИТЕРАТУРЫ

1. Bistatic radar: emerging technology / ed. by M. Cher-niakov. Chichester: John Wiley & Sons, 2008. 406 p.

2. Сергиенко А. Б. Цифровая обработка сигналов: учеб. пособие. 3-е изд. СПб.: БХВ-Петербург, 2011. 768 с.

3. Altera. May 2013. Radar Processing: FPGA or GPUs? WP-01197-2.0 White paper. URL: http://www.altera.com /literature/wp/wp-01197-radar-fpga-or-gpu.pdf (дата обращения: 03.09.2015).

A. V. Barkhatov, A. S. Kozlov Saint Petersburg Electrotechnical University "LETI"

Radar amplitude-range-doppler surface fast calculation on graphics processing units

The paper describes the algorithm of fast calculation of the amplitude-range-Doppler surface. The features of the software implementation of the algorithm are drawn. The nuances of parallel computing of the surface on graphics processing units of the gaming graphics card are shown. The result, the sixteen surfaces parallel computing time for multichannel radar, is presented.

Radar, digital processing, ambiguity function, cross-correlation, Doppler shift, delay, fast Fourier transformation, graphics processing units, parallel computing

Статья поступила в редакцию 7 сентября 2015 г.

УДК 621.396.9

С. В. Штанько, Д. А. Лесняк Военно-космическая академия им. А. Ф. Можайского

Алгоритмы защищенного информационного обмена в радиоканалах космической навигационной системы

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

Космическая навигационная система, межспутниковые каналы, алгоритм Диффи-Хеллмана, криптографическая защита

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

На современном этапе развития радиотехнических систем особое место занимают космические навигационные системы (КНС).

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

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

аппаратура потребителя (НАП)", что накладывает существенные ограничения на процедуры аутентификации и обмена ключевой информацией.

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

© Штанько С. В., Лесняк Д. А., 2015

47