МЕТОД РАСЧёТА СЕМЕЙСТВА ФУНКЦИЙ БЕССЕЛЯ С ПОМОЩЬЮ БПФ ПРИ РЕАЛИЗАЦИИ ТЕОРЕТИЧЕСКИХ ОСНОВ АКТИВНО-ПАССИВНОЙ НИЗКОЧАСТОТНОЙ ТОМОГРАФИИ Текст научной статьи по специальности «Физика»

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

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

Исследования СОК позволили выявить целый ряд ее преимуществ по сравнению с ПСС, которые приведены в [4, 6]. Основным недостатком СОК является сложность реализации немодульных операций, таких как перевод чисел в СОК и обратно и других. Однако базовые операции цифровой фильтрации являются модульными и могут быть выполнены за один модульный такт (за время выборки результата операции из памяти ЦПОС), которое в настоящее время составляет от 0,1 до 10 нс. Алгоритм цифровой нерекурсивной фильтрации идеальным образом согласуется с организацией вычислений в СОК [4, 7].

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

Проведенные исследования показали, что кроме значительного повышения быстродействия реализации алгоритма цифровой фильтрации применение СОК позволило повысить прямоугольность амплитудно-частотной характеристики цифрового фильтра на 13 % и точность вычислений выходных отсчетов на порядок [7].

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

Литература

1. Солонина А.И., Улахович Д.А., Яковлев Л.А. Алгоритмы и процессоры обработки сигналов. СПб., 2002.

2. Цифровая обработка сигналов: Учебн. пособие для вузов / Л.М. Гольденберг, Б.Д. Матюшкин, М.Н. Поляк: 2-е изд. М., 1990.

3. Червяков Н.И., Велигоша А.В., Калмыков И.А. Отказоустойчивый цифровой фильтр обработки комплексных данных. Ставрополь, 1995. С. 82-83.

4. Коляда А.П., Пак И.Т. Модулярные структуры конвейерной обработки цифровой информации. Минск, 1992.

5. Нейрокомпьютеры в остаточных классах. Кн. 11 Н.И. Червяков, П.А. Сахнюк, А.В. Шапошников): Учеб. пособие для вузов. М., 2003.

6. Eyre J., Bier J. The Evolution of DSP Processor//IEEE Signal Processing magazine. 2000. March.

7. Червяков Н.И., Велигоша А.В., Калмыков И.А., Иванов П.Е. Цифровые фильтры в системах остаточных классов // Радиоэлектроника. 1995. Т. 38. № 8. С. 11-20.

21 февраля 2006 г.

УДК 621.391.26

МЕТОД РАСЧЁТА СЕМЕЙСТВА ФУНКЦИЙ БЕССЕЛЯ С ПОМОЩЬЮ БПФ ПРИ РЕАЛИЗАЦИИ ТЕОРЕТИЧЕСКИХ ОСНОВ АКТИВНО - ПАССИВНОЙ НИЗКОЧАСТОТНОЙ ТОМОГРАФИИ

© 2006 г. П.А. Стародубцев, А.Н. Иванов, А.А. Жевагин

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

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

Поэтому в 70-80 гг. прошлого столетия специалистами Акустического института имени академика

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

Но как показала практика, основными недостатками параметрических просветных методов, ограничивающими их широкое применение для обнаружения ПО, являются:

- низкая эффективность параметрического преобразования при нелинейном взаимодействии исходных высокочастотных волн в низкочастотные комбинационные;

- интенсивность волн комбинационных частот зависит от нелинейных характеристик среды (солености, температуры, насыщенности микроорганизмами или газовыми пузырьками) и изменяется в широких пределах от 4 -10-7 до 19 -10-5 относительно интенсивности высокочастотных волн накачки;

- накапливающееся с расстоянием комбинационное взаимодействие сферических волн, рассеянных отдельными рассеивателями на частоте ПС с квазиплоской, инфранизкочастотной (ИНЧ) волной шумоиз-лучения ПО наблюдается достаточно слабо.

В классических параметрических системах данное взаимодействие составляет 0,01... 0,1 %. В системах с применением дополнительных мер и устройств (ультразвуковая пузырьковая область) эффективность преобразования энергии данных волн может составлять 0,1 %.

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

- низкочастотном просветном методе широкомасштабного интенсивного нелинейного взаимодействия и параметрического преобразования ПС возмущенной областью, созданной движущемуся ПО во всех направлениях их совместного распространения;

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

- гидроакустической системе как широкомасштабном параметрическом излучателе-приемнике с низкочастотной подсветкой (накачкой) контролируемой трассы, являющейся одновременно и объемной гидроакустической базой.

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

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

2. Нелинейное взаимодействие несущей низкочастотного ПС с наиболее интенсивными ИНЧ составляющими возмущенной области, созданной движущимся ПО, и параметрического преобразования его в комбинационные волны суммарной и разностной частоты.

Механизм генерации данных физических явлений определяется как монопольными акустическими колебаниями элементарных рассеивателей водной среды, так и перемещением ее деформаций как устойчивой совокупности смещений частиц пограничного слоя воды вокруг корпуса ПО и его окружения. При этом перерассеяние и изменение параметров ПС происходит во множестве точек и слоев перемещающейся стратифицированной водной среды во всем объеме и во всех направлениях за счет возникновения внутренних корабельных волн. На глубинах меньших 100 м отношение амплитуд внутренних волн к глубине (параметр нелинейности) в этом случае достигает значения 0,1 и более. Возмущения гидрологических полей, в том числе скорости звука, таким волновым движением будут качественно отличаться от аналогичных возмущений, создаваемых акустическими с бесконечно малыми амплитудами, т. е. линейными волнами от шума корпуса ПО. При этом различные части корпуса ПО регенерируют гармонические и квазигармонические полнопериодные, полупериодные и четверть периодные ИНЧ волны (О) в вертикальной и горизонтальной плоскости при неустойчивом движении ПО. При распространении такой суммарной внутренней волны с конечными амплитудами происходит постепенное изменение формы волны вследствие разности в скоростях движения различных участков его профиля. Периодические волны в процессе распространения эволюционируют в возмущения с «бо-рообразным» (подобным профилю ударных волн) профилем.

Увеличение крутизны внутренней волны соответствует нарастанию высокочастотных гармоник, т. е. обогащению ее спектра высокочастотными составляющими волнового поля, которые вызывают интенсивное изменения спектра ПС. Спектр такого ПС в общем случае содержит бесконечное число гармонических составляющих, частоты которых равны ю 0 ± кО, а интенсивность этих составляющих пропорциональны значениям функции Бесселя 3 к (тк).

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

щенной области ПЛ низкочастотный ПС представляет собой фазомодулированную волну в пределах одной лучевой трубки, суммарное давление в которой на приемнике с учетом гидрологических условий по трассе распространения ПС определяется следующей математической зависимостью:

Ру (Ь-, е г ) = [Р0(0« сп х

от выражения (1), функция Бесселя 3п(х), формула которой для целого и нецелого значения порядка п имеет следующий вид [1]:

Í „2 Ak

Jn (x) =

kfo k!Г(n + k +1)

(2)

X( E Jk (mk )C0S Ю 0

k=1

tt - -

c0

+ kQ, (1 - cose г )ttmk)] ± P±E, (1)

где ti - суммарное время распространения ПС по разным акустическим лучам в пределах / лучевой трубки; асп - коэффициент пространственного затухания на частоте ПС; — - время фазового модуляци-

с 0

онного (параметрического) взаимодействия низкочастотного ПС с ИНЧ гидрофизическими и другими полями возмущенной области, созданной движущемся ПО, в пределах /-ой лучевой трубки; к - любое вещественное число от 1 до / ; тк - индекс фазовой модуляции,

где Г (п) = | е Чп 1С( - гамма-функция, которая при

0

целочисленном значении аргумента п может быть интерполирована с помощью выражения Г(п+1)=п!. Данная функция является решением дифференциального уравнения Бесселя для целого и нецелого значения п:

d 2 y

dy

-Г + x

dx2 dx

dx +(x2 + n2 )y = 0.

В литературе [1] встречается другое выражение для расчёта функции Бесселя

1 2п

Jn (x) = — J exp [jx sin 9 - jn9] d9 =

2n

mk =

(Y-1+cos e,. )P0 (t)

2po(co)3P (Q,,Ьг)

sin

n L- (1 - cos e,) Q, ,

П -i- (1 - cos e, ) Q, ,

1 2n

= — J exp [jx sin e]exp [- jne] d e, 2п 0

(3)

Аю Q,

которое по существу представляет собой преобразование Фурье от функции ехр[/х^ше], где х - постоянная величина, е - переменная эквивалентная времени, п - переменная эквивалентная частоте.

2П

Выражение (3) с помощью замены Се ^ — и

3 к (тк) - функция Бесселя к индекса от аргумента т к , показывающая зависимость коэффициента модуляции от амплитуды модулирующего сигнала; ю - е ^ к— может быть преобразовано в формулу ДПФ

N

2п

несущая частота ПС; Аю - девиация частоты ПС; О / -ИНЧ колебания возмущенной области ПО с частотами 0О 2,...,О/ соответственно; Р0(/) - исходное давление ПС на излучателе; Ь/ - проекция плоскости сечения области взаимодействия (в том числе и параметрического с появлением комбинационных волн разностной и суммарной частоты) ПС с ИНЧ составляющими возмущенной области на линию, соединяющую излучатель и приемник; Р(0 /, Ь/) - суммарное давление /-й гармоники ИНЧ в возмущенной области, созданной движущемся ПО; р 0 - плотность водной среды; е / - угол между взаимодействующими ПС и гармониками ИНЧ возмущенной области; у -параметр нелинейности среды (морская вода - 3,5 ед., возмущенная область - 10 .15 ед.); ± Р±у - давление

волн комбинационных частот суммарной и разностной частот.

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

N

следующего вида:

1 N -1

Jn(x) = N £ exP

N k=0

. , 2п jx sinl — k

1 N

exp

2n , - j—nk

N

1 N -1

— E exp

N k=0

2n , jx sinl — k

1 N

WN

(4)

где Wn = exp

2n , - j—nk

N

- комплексная величина

эквивалентная фазовому коэффициенту в формуле ДПФ; к = 0, ..., N - 1 - номер выборки по е.

В ряде задач, где используется функция Бесселя 3п(х) [2], требуется выполнять обработку информации в реальном масштабе времени. Для расчета функции с помощью выражения (4) в данной статье предлагается использовать один из известных алгоритмов БПФ, например, алгоритм БПФ по основанию 2 (алгоритм Кули - Тьюки). В этом случае, для одного из N отсчётов переменной х, будут рассчитываться сразу N значений функций 3п(х) порядков п = 0, ..., N - 1.

Результаты подобных расчетов в системе инженерных и научных расчётов МАТЬАВ 6.5 отображены на рис. 1.

0

Jn(x)

-0,5

10 12

14

б)

0,4 0,3 0,2 0,1 0 -0,1 -0,2 -0,3

AJn(x)

0 2

10 12

14

в)

AJn(x)x10-1

0,8 0,6 0,4 0,2 0 -0,2 -0,4 -0,6 -0,8

В вычислительном эксперименте выполнялись расчёты функций Бесселя для различных порядков n = 0, ..., 15 (с шагом 1) для области определения х = 0, ..., 16 (с шагом 0,01). В самом начале выполнялись расчёты функций Бесселя по традиционным соотношениям (1), результат вычислений с помощью функции MATLAB besselj показан на рис. 1а. Затем выполнялся расчёт с помощью алгоритма БПФ (функция MATLAB fft), рис. 16. На указанных рисунках показаны графики функции Бесселя только для порядка n = 0, ..., 7. После этого результаты вычислений сравнивались между собой. На рис. 1 в и г показаны графики абсолютной ошибки вычислений для различных диапазонов х= 0, ..., 16 и х=0, ..., 4, которая рассчитывалась как разница AJn(x) между результатами двух различных методов.

Как видно из графиков на рис. 1 в и г, разница растёт быстрее с увеличением порядка n, а при n=7 наблюдается самый быстрый её рост. Но даже для больших порядков на рис. 1г ошибка вычислений не превышает 10 - 3 для х = 0, ..., 4. Такую особенность данного метода вычисления можно объяснить двумя причинами: периодичностью спектра при вычислении ДПФ и явлением наложения (aliasing), наблюдаемую при неправильной дискретизации непрерывных величин [2].

Известно, что результат дискретного преобразования Фурье имеет периодическую и симметричную структуру, поэтому из него надо выбрать только часть составляющих для n=0, ..., N/2 - 1. В выражении (2)

• I 2п

гармоническая функция exp

jx sin

N

зависит от

двух переменных: х, которая играет роль частоты, sin I — к I, которое представляет собой периодически

0 0,5 1,0 1,5 2,0 2,5 3,0 3,5 х

г)

Рис. 1. Результаты расчётов в системе ЫАТЬАВ: а - по выражению (1), б - расчёт с помощью БПФ, в и г - ошибка вычисления для х = 0, и х = 0, ..., N/4

N

изменяющееся в диапазоне от -1 до 1 дискретное время. Поэтому при росте переменной х выше N/2 происходит наложение периодических составляющих преобразования Фурье. Для иллюстрации приведём трёхмерные графики функции Бесселя на рис. 2, рассчитанные в ЫАТЬАВ рассматриваемыми выше методами для значений порядка п=0, ..., 63 (с шагом 1) и независимой величины х=0, ..., 31 (с шагом 0,1).

Из рис. 2в видно, что различие между результатами расчётов равно нулю, если порядок п=0, ..., 31, а независимая переменная х = 0, ..., 31. На основании вышесказанного следует, что при вычислении функций Бесселя для порядков п = 0, ..., N/2 следует рассчитывать с помощью алгоритма БПФ сразу N значений, но учитывать как результат только первые N/2 значения, потому что остальные составляющие будут периодически повторяться. Кроме этого, следует помнить, что эти вычисления более точны методом БПФ только при х < N/2 (в данном случае играет роль частоты Найквиста). При х > N/2 абсолютное значение ошибки будет резко увеличиваться из-за эффекта наложения, наблюдаемого в ДПФ. Но при проведении расчетов для низкочастотной акустической томографии данных точностных характеристик вполне достаточно.

0

2

4

6

8

x

4

6

8

x

7

6

Примечание.* - следует учесть, что предлагаемый метод расчёта (с помощью БПФ) требует в два раза больше значений по п=0, ..., 2Г - 1, т.е. г = 3, ..., 11 вместо г = 2, .,10.

Таблица

Время расчёта семейства функций Бесселя для п = 0, ..., 2Г- 1*

Метод расчёта r

2 3 4 5 б 7 8 9 10

Время расчёта, с

Традиционный 0,12 0,151 0,29 0,521 0,8б1 1,4б2 3,б75 8,122 25,397

С помощью БПФ 0,37 0,821 1,783 5,047 10,7бб 21,7б1 43,342 8б,405 172,39

в)

Рис. 2. К описанию особенностей вычисления функций Бесселя: а - традиционным методом; б - с помощью БПФ; в - разница между ними

Для сравнительной оценки быстродействия предлагаемого метода расчёта в системе MATLAB 6.5 проводился вычислительный эксперимент. С помощью

Морской государственный университет им. адмирала Г.И. Невельского, г. Владивосток

функций программного таймера tic и toc, описанных в [3], определялось время, затраченное на вычисления. Полученные при этом результаты для традиционного [1] и предлагаемого методов расчёта приведены в таблице и на рис. 3. Они показывают, что предлагаемый метод имеет в 3...8 меньшее быстродействие, чем расчёт функции Бесселя традиционным способом по выражению (2).

Время, с

160 120 80 40

0 2 3 4 5 6 7 8 9 2Г

Рис. 3. К сравнительному анализу методов расчёта: 1 - расчёт традиционным способом; 2 - расчёт с помощью БПФ

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

Литература

1. Математика: Большой энциклопедический словарь/ Гл. ред. Ю.В. Прохоров: 3-е изд. М., 1998.

2. Рабинер Л., Гоулд Б. Теория и применение цифровой обработки сигналов/ Пер. с англ. А.Л. Зайцева, Э.Г. На-заренко, Н.Н. Тетёкина. М., 1978.

3. Потёмкин В.Г. Система инженерных и научных расчётов MATLAB 5.x: В 2 т. М., 1999.

31 мая 2006 г.