CC BY
76Электронный журнал «Техническая акустика» http://www .ejta.org
2012, 5
А. А. Хамухин
Национальный исследовательский Томский политехнический университет, Россия, 634050 Томск, пр-т Ленина, 30 e-mail: aaxtpu@tpu.ru
О параллельном вычислении непрерывного вейвлет-преобразования в задаче обнаружения узкополосных шумовых гидроакустических сигналов на основе интегрального вейвлет-спектра
Получена 24.08.2012, опубликована 25.09.2012
Установлено, что в задаче обнаружения узкополосных шумовых гидроакустических сигналов на основе интегрального вейвлет-спектра значительное время занимает вычисление непрерывного вейвлет-преобразования, что снижает эффективность обнаружения. Предложен алгоритм параллельного вычисления непрерывного вейвлет-преобразования оцифрованных отсчетов шумовых сигналов, позволяющий сокращать время вычисления интегрального вейвлет-спектра. Показано, что сокращение времени достигается за счет предварительного расчета матриц базисных вейвлет-коэффициентов и сохранения их в запоминающем устройстве, используемом при дальнейшей онлайн-обработке сигналов. Представлены результаты компьютерного моделирования обнаружения узкополосных шумовых гидроакустических сигналов с помощью предложенного алгоритма.
Ключевые слова: узкополосный шумовой гидроакустический сигнал, интегральный вейвлет-спектр, непрерывное вейвлет-преобразование, параллельные вычисления
ВВЕДЕНИЕ
Задача обнаружения узкополосных шумовых гидроакустических сигналов, известная также как задача обнаружения дискретных составляющих шумового спектра движущегося объекта, в настоящее время решается различными способами [1-4]. Известно, что шум, излучаемый движущимся в воде объектом, является смесью шумов двух типов: широкополосный шум с непрерывным спектром и узкополосный шум с дискретным спектром, состоящий из синусоидальных составляющих на отдельных частотах (дискретные составляющие). Характеристики дискретных составляющих (ДС), обусловленные работой вспомогательных механизмов объекта, обычно стабильны и не зависят от скорости движения объекта. Частоты и амплитуды ДС, вызванные энергетической установкой и гребными винтами объекта, могут изменяться в зависимости от скорости и глубины погружения объекта. ДС от некоторых
источников образуют звукоряды (с кратными частотами), которые являются высокостабильными и имеют узкие полосы частот [4].
Эти факты усложняют задачу обнаружения реальных объектов и обуславливают проведение новых исследований и применение новых методов.
Так, например, метод непрерывного вейвлет-преобразования (НВП), появившийся относительно недавно [5], находит применение и в задачах акустики [6-9]. НВП является мощным инструментом анализа нестационарных периодических сигналов, позволяющим, в отличие от преобразования Фурье, получать частотно-временную характеристику сигнала. На основе этой характеристики можно создавать и устройства для обнаружения узкополосных шумовых гидроакустических сигналов.
Так, в работе [4] описывается подобное устройство, в котором после измерения и оцифровки исследуемого сигнала вычисляется его интегральный вейвлет-спектр, по которому пороговое устройство выдает сигнал обнаружения дискретной составляющей. Основная вычислительноёмкая операция этого устройства -непрерывное вейвлет-преобразование (ВП), время расчета которого увеличивается экспоненциально с ростом количества отсчетов в выборке анализируемого сигнала. Для ускорения вычисления вейвлет-преобразования разработано множество методов: ВП на основе быстрого преобразования Фурье, дискретное ВП, диадное ВП и др. [10]. Но все эти методы связаны с частичной потерей информации при вычислении, что может снизить эффективность обнаружения узкополосных шумовых гидроакустических сигналов. Поэтому для данной задачи важно сохранить классическое непрерывное вейвлет-преобразование и сокращать время его расчета за счет распараллеливания вычисления НВП, что и является целью настоящей работы.
1. ДИСКРЕТИЗИРОВАННОЕ НЕПРЕРЫВНОЕ ВЕЙВЛЕТ-ПРЕОБРАЗОВАНИЕ
Классическое непрерывное вейвлет-преобразование произвольного непрерывного сигнала S(t) описывается следующим уравнением [10]:
(материнский вейвлет); Ь — сдвиг по времени; a — масштаб по времени (обратно пропорциональный частоте); W(a, Ь) — искомое вейвлет-преобразование.
При цифровой обработке непрерывный сигнал S(t) заменяется последовательностью его дискретных отсчетов, измеренных с частотой, удовлетворяющей теореме Котельникова-Найквиста. Поэтому на практике реально вычисляется дискретизированное непрерывное вейвлет-преобразование, за которым далее сохраним термин НВП и обозначение W(a, Ь), чтобы отличать его от дискретного ВП. В соответствии с этим интеграл в уравнении (1) аппроксимируем конечной суммой и получим:
(1)
где
некоторая функция с определенными свойствами
Ж (а, Ь) к Ж (а}., Ьк );
(2)
где г, / к — индексы по времени t, по масштабу а, по сдвигу по времени Ь; N — количество шагов по времени; М — количество масштабов а/; Дti — шаг по времени; S(tг■) — оцифрованный сигнал.
В работах [11, 12] описано, как формулу (2) можно реализовать с помощью ячеек однородной вычислительной структуры. При этом достигается максимально возможное синхронное распараллеливание операций. Это распараллеливание не на 4 или 16 ядер современных микропроцессоров, а на тысячи и более однотипных параллельно работающих вычислительных ячеек, которые могут быть реализованы в виде схемы, приведенной в работе [13].
Однако создание однородных вычислительных структур является самостоятельной технологической проблемой, а алгоритм параллельного вычисления НВП можно применять уже сейчас на имеющихся вычислительных средствах. Это не позволит достигнуть такого кардинального ускорения расчета НВП, как на однородных вычислительных структурах, но все же это тоже ускорение. Для этого алгоритм, представленный в работах [11-13] необходимо модифицировать.
1.1. Расчетные формулы вычисления непрерывного вейвлет-преобразования
Распараллеливание вычислений в настоящее время является основным резервом повышения производительности вычислительных средств. Большинство задач имеет свои особенности, требующие разработки специальных алгоритмов распараллеливания. Есть эти особенности и при вычислении вейвлет-преобразования. Наиболее полный обзор параллельно-рекурсивных алгоритмов вычисления ВП приведен в работе [14]. Однако они в основном предназначены для распараллеливания дискретного вейвлет-преобразования, широко применяемого для обработки видео- и аудиоинформации.
Для распараллеливания непрерывного ВП возьмем за основу рекуррентную (по индексу г) формулу, предложенную для параллельного вычисления НВП в ячейках однородных вычислительных структур [11]:
входной сигнал S(t), (/=0-^; Рг/л — коэффициенты, полученные из материнского вейвлета.
Для создания параллельного алгоритма вычисления НВП предлагается задать коэффициенты Рг/к в несколько ином виде, чем они представлены в работах [11, 12]:
Ж1 , Ьк) = Ж1 а, Ьк) + 5 (Ш). ,
(3)
где Ж0(а/, Ьк)=0 — начальное условие; S(iДf) — оцифрованный с равномерным шагом Дt
р, 1,к =
где ¥
У
(і - к )Аї
а.
(4)
(і - к )Аї
— некоторая функция с определенными свойствами, полученная из
материнского вейвлета.
Если сдвиг по времени, по которому вычисляется вейвлет преобразование, принять равным шагу оцифровки входного сигнала Д^, то тогда Ьк=кА1. В этом случае, учитывая, что 1г=1А1, в уравнении (4) аргумент функции ¥ можно преобразовать, и мы получим:
У
(і - к )Аґ
а
= У
(5)
Тогда формулу (3), взятую в качестве исходной, можно записать в виде:
Ж1+1 (а і, Ьк ) = Ж1 (а,, Ьк ) + 5(іМ)
У
(* А ^
1г - Ьк
а,
V 1 J
(6)
Раскрывая рекурсию по индексу г в выражении (6), мы получим точное совпадение с формулой (2) для дискретизированного непрерывного вейвлет-преобразования при условии оцифровки сигнала с равномерным шагом по времени. Таким образом, мы доказали, что выражение (3) с модифицированным видом вейвлет-коэффициентов (4) можно использовать для построения параллельно-рекурсивного алгоритма вычисления дискретизированного непрерывного вейвлет преобразования.
1.2. Алгоритм параллельного вычисления непрерывного вейвлет-
преобразования
Распараллеливание вычисления формулы (3) можно выполнять и по масштабам а/, и по сдвигам по времени Ьк. На схему устройства с полным распраллеливанием по а/ и Ьк получен патент РФ [13]. Однако такой параллельный алгоритм требует большого количества процессоров в вычислительном устройстве. Поэтому представляется актуальным преобразовать расчетные формулы к такому виду, чтобы можно было применять частичное распараллеливание в соответствии с имеющимся в распоряжении количеством процессоров.
Преобразуем формулу (3) для распараллеливания по масштабам а/.
Ж^Ц, Ьк) = Жг (аі, Ьк) + 5 (іМ). РіХк;
Жг+1(а2, Ьк) = Жі (а2, Ьк) + 5 (іМ). р
(7)
Жм(ам, Ьк) = Жі (ам, Ьк) + 5 (іМ) • Рг
, м , к
Все уравнения выражения (7) независимы друг от друга и их можно вычислять параллельно.
Формула (3) для распараллеливания по сдвигам по времени Ьк будет иметь вид: Жм(а], Ь) = Жг (а/, Ьк ) + 5 (гМ) • Р,; л;
Жг+1 (а/, Ь2) = Ж- (а/, Ьк ) + 5 (гМ) • Рии 2; (8)
Жм(а], ЬN ) = Жг (а/, Ьк ) + 5 (гМ) • Ри,„.
Выражение (8) также состоит из независящих друг от друга уравнений, которые можно вычислять параллельно.
Коэффициенты Рг/к в выражении (3) и в уравнениях (7) и (8) не зависят от входного сигнала, поэтому их можно рассчитать заранее для всего диапазона индексов г, / к. Для этого только необходимо выбрать вид функции ¥ и иметь достаточный объем памяти для хранения этих коэффициентов, что в настоящее время не является проблемой. Например, если взять выборку из N=2000 отсчетов и М=600 масштабов (как в работе [4]), то объем требуемой памяти для хранения массива 32-разрядных коэффициентов Рг/к составит около 8 Гбайт. Такой подход предварительного расчета матриц базисных вейвлет-коэффициентов Рг/к дает основной вклад в сокращение времени расчета НВП.
Формулы для параллельного вычисления НВП (7) и (8) имеют рекуррентный вид и могут вычисляться итерационно. Это значит, что рассчитанные на предыдущем шаге значения Ж(а/, Ьк) сохраняются в запоминающем устройстве и используются при поступлении следующего отсчета входных данных. Это дает дополнительный вклад в сокращение времени расчета НВП.
Таким образом, весь алгоритм вычисления дискретизированного непрерывного вейвлет-преобразования разделяется на два этапа.
1. Подготовительный этап. Для успешной реализации этого этапа необходимо иметь некоторую предварительную информацию об анализируемом сигнале. С помощью этой информации выбираются: шаг оцифровки анализируемых данных Д(, количество и значения масштабов по времени (М и а/), вид функции ¥ (материнский вейвлет), ограничение на максимальное количество отсчетов анализируемого сигнала (Щ. Затем по формуле (4) вычисляются и сохраняются в запоминающем устройстве все базисные коэффициенты Ру,к , которые потребуются для реализации второго этапа вычисления НВП.
2. Основной этап. По формулам (7) или (8) параллельно вычисляется текущее значение НВП для текущего момента времени (гД). При поступлении следующего отсчета входных данных в момент времени (г'+1)Д^ к уже рассчитанному массиву значений НВП добавляется только один слой значений, рассчитываемых по формулам (7) или (8). При этом используется подготовленный на первом этапе массив коэффициентов Рг/к . Таким образом, реализуется итерационный процесс уточнения НВП, что позволяет уменьшить количество вычислений и сократить время расчета.
Следует отметить, что большинство вычислительных операций переносится на подготовительный этап, результатом которого является достаточно большой массив чисел (коэффициенты Р/к). На втором этапе вычислительные затраты относительно невелики, но требуется достаточно много памяти для хранения подготовленного на первом этапе массива чисел.
2. КОМПЬЮТЕРНОЕ МОДЕЛИРОВАНИЕ ОБНАРУЖЕНИЯ УЗКОПОЛОСНЫХ ШУМОВЫХ ГИДРОАКУСТИЧЕСКИХ СИГНАЛОВ С ПОМОЩЬЮ ПРЕДЛОЖЕННОГО АЛГОРИТМА
Задача обнаружения узкополосных шумовых гидроакустических сигналов на основе интегрального вейвлет-спектра подробно изложена в работе [4].
Интегральный вейвлет-спектр Щ(а), называемый еще скалограммой, рассчитывается на основе вычисленного непрерывного вейвлет-преобразования и отражает вклад различных масштабов (частот) в полную энергию входного сигнала. Наличие максимумов в скалограмме позволяет выявлять частоты, которые превышают общий уровень, и таким образом обнаруживать дискретные составляющие шумовых гидроакустических сигналов.
Для непрерывного сигнала интегральный вейвлет-спектр имеет вид [4]:
где Т — полный интервал времени наблюдения входного сигнала.
Для дискретизированного входного сигнала интеграл в выражении (9) аппроксимируется суммой:
Формулу (10) необходимо вычислять каждый раз после вычисления НВП по любой из приведенных выше формул.
Компьютерная модель представленного здесь алгоритма создана в среде Excel с программированием на языке VBA. Для тестирования алгоритма был выбран входной сигнал, содержащий две ДС на фоне аддитивной помехи. Было выбрано количество отсчетов входного сигнала N=256, количество масштабов M=16, шаг Д/=1. Реальная задача обнаружения узкополосных шумовых гидроакустических сигналов на основе интегрального вейвлет-спектра решается при несколько больших количествах ДС, масштабов и отсчетов, но для проверочного тестирования работоспособности алгоритма мы сочли эти значения достаточными.
Тестирование проводилось на двух материнских вейвлетах: вейвлете Морле и вейвлете «мексиканская шляпа» (MHAT) [10]. Их вид представлен на рис. 1.
0
(9)
(10)
Рис. 1. Графический вид вейвлетов: Морле (слева) и МНАТ (а=4, Ь=50)
Входной сигнал 5(і) моделировался в виде суммы двух гармоник (с частотами 1 и 4 Гц) с помехой в виде модели дискретного белого шума. Вид анализируемых входных сигналов с различным уровнем помехи представлен на рис. 2.
Рис. 2. Входной сигнал с различным уровнем помехи (а - без помехи, Ь - 12 дБ, с - 16 дБ)
Вид получаемого интегрального вейвлет-спектра зависел от вида применяемого материнского вейвлета, его параметров и уровня помех. Некоторые виды расчетного интегрального вейвлет-спектра приведены на рис. 3.
#0)80
Масштаб (а)
ТУ (а)80
1 3 5 7 9 11 13 15
Масштаб (а)
Ф(а)80
1 3 5 7 9 11 13 15
Масштаб (а)
с)
Рис. 3. Расчетный интегральный вейвлет-спектр для входных сигналов с различным уровнем помехи (а - без помехи, Ь - 12 дБ, с - 16 дБ)
Как видно из графиков, расчетный интегральный вейвлет-спектр содержит два максимума, соответствующих двум гармоникам во входном сигнале. С увеличением уровня помехи выраженность максимумов несколько уменьшается, и при уровне помехи в 16 дБ второго максимума уже практически нет, т. е. низкочастотная
дискретная составляющая спектра входного сигнала не обнаруживается. Такой
результат характерен в целом для метода обнаружения дискретных составляющих на основе интегрального вейвлет-спектра и не является следствием параллельного
вычисления непрерывного вейвлет-преобразования по предложенному в данной работе алгоритму.
Для приведенных выше исходных данных время расчета массива Pij,k составляло около 30 с (подготовительный этап предлагаемого алгоритма), а время расчета НВП составляло около 0,5 с на 1 масштаб (основной этап предлагаемого алгоритма). Для 16 масштабов время расчета составляло около 8 секунд. Таким образом, в период самого процесса обнаружения время расчета сократилось примерно в 2 раза за счет вынесения расчетов матриц базисных коэффициентов в подготовительный этап. Расчеты производились на микропроцессоре Intel Core i7 3,40 GHz под управлением Windows 7(64). При увеличении значений N и M выигрыш по времени расчета основного этапа алгоритма будет больше.
ЗАКЛЮЧЕНИЕ
В задаче обнаружения узкополосных шумовых гидроакустических сигналов на основе интегрального вейвлет-спектра значительное время занимает вычисление непрерывного вейвлет-преобразования, что снижает эффективность обнаружения. Предложенный алгоритм параллельного вычисления непрерывного вейвлет-преобразования оцифрованных отсчетов шумовых сигналов, позволяет сокращать время вычисления интегрального вейвлет-спектра. Сокращение времени достигается за счет предварительного расчета матриц базисных вейвлет-коэффициентов и сохранения их в запоминающем устройстве, используемом при дальнейшей онлайн-обработке сигналов. Результаты компьютерного моделирования обнаружения узкополосных шумовых гидроакустических сигналов подтвердили работоспособность предложенного алгоритма. Предложенный алгоритм можно реализовать как программно на универсальных микропроцессорах, так и в виде специализированных вычислительных устройств. Аппаратная схема реализации предложенного алгоритма имеет патент РФ № 2437147.
ЛИТЕРАТУРА
1. Калью В. А. Автоматический мониторинг подводного шума судов с помощью искусственной нейронной сети. [Электронный ресурс] // Российский фонд фундаментальных исследований, грант № 11-08-00402. Режим доступа: http://www.rfbr.ru/rffi/ru/contest/n_549/o_41432, свободный.
2. Калью В. А., Неворотин В. Ю., Правдин А. А. Способ обнаружения сигнала источника, порождающего дискретную составляющую в спектре суммарного сигнала нескольких источников: пат. 2393490 Рос. Федерация. № 2009103585/28; заявл. 03.02.2009; опубл.: 27.06.2010, Бюл. № 18. - 11 с.
3. Калью В. А., Неворотин В. Ю., Правдин А. А. Повышение точности локализации дискретных составляющих спектра шума движущегося транспортного средства. [Электронный ресурс] // Техническая акустика. - Электрон. журн. - 2012. - 3.
Режим доступа: http://www.ejta.org, свободный.
4. Сапрыкин В. А., Малый В. В., Шаталов Г. В. Устройство обнаружения узкополосных шумовых гидроакустических сигналов на основе вычисления интегрального вейвлет-спектра: пат. 2367970 Рос. Федерация. № 2007145474/28; заявл. 28.11.2007; опубл. 20.09.2009, Бюл. № 26. - 27 с.
5. Daubechies I. The wavelet transform, time-frequency localization and signal analysis // IEEE Trans. Inform. Theory.- 1990, №5, p. 961-1005.
6. Седдики А., Джебари А., Руваен Ж. М. Использование стационарного вейвлет-преобразования в БЧХ-кодировании изображений. [Электронный ресурс] // Техническая акустика. - Электрон. журн. - 2003. - 1. Режим доступа: http://www.ejta.org, свободный.
7. Джеббоури М., Джебури Д., Наом Р. Применение Фурье- и вэйвлет-преобразований для восстановления «зашумленного» изображения. [Электронный ресурс] // Техническая акустика. - Электрон. журн. - 2003. - 2. Режим доступа: http://www.ejta.org, свободный.
8. Бограчев К. М. Сравнение эффективности Фурье- и вейвлет-декомпозиции в пассивной акустической термотомографии // Акуст. журн., 2005, т. 51, № 3, с. 239245.
9. Эль Аллами М., Рхимини Х., Насим А., Сидки М. Применение вейвлет-преобразования к определению характеристик волн Лэмба. [Электронный ресурс] // Техническая акустика. - Электрон. журн. - 2010. - 8. Режим доступа: http://www.ejta.org, свободный.
10. Яковлев А. Н. Введение в вейвлет-преобразования. - Новосибирск: Изд-во НГТУ, 2003. - 104 с.
11. Хамухин А. А. Применение ячеек однородной структуры для вычисления непрерывного вейвлет-преобразования // Известия Томского политехнического университета. - 2010. - Т. 317. - № 5. - С. 149-153.
12. Хамухин А. А. Математическая модель ячейки однородной структуры для вычисления непрерывного вейвлет-преобразования // Проблемы информатики. -2011. - № 5. - С. 87-93.
13. Хамухин А. А. Устройство для вычисления дискретизированного непрерывного вейвлет-преобразования: пат. 2437147 Рос. Федерация. № 2010127068/08; заявл. 01.07.2010; опубл. 20.12.2011, Бюл. № 35.- 9 с.
14. Нго Хыу Фук. Параллельно-рекурсивные методы выполнения вейвлет-преобразования в задачах обработки дискретных сигналов // Электронная библиотека диссертаций. 2005. URL: http://www.dissercat.com/content/parallelno-rekursivnye-metody-vypolneniya-veivlet-preobrazovaniya-v-zadachakh-obrabotki-disk (дата обращения: 11.05.2012).
