Эффективность вейвлет фильтрации сигнала на GPGPU (стр. 2 из 6)
Рисунок 3 - Трехуровневый банк фильтров
На каждом уровне вышеприведённой диаграммы сигнал раскладывается на низкие и высокие частоты. В силу двукратного прореживания длина сигнала должна быть кратна 2n, где n — число уровней разложения.
Например, для сигнала из 32 отсчётов с частотным диапазоном от 0 до fn трёхуровневое разложение даст 4 выходных сигнала в разных масштабах:
Таблица 1– Трехуровневое разложение сигнала размером в 32 отсчета
| Уровень | Частоты | Длина сигнала |
| 3 | 0 … fn / 8 | 4 |
| fn / 8 … fn / 4 | 4 | |
| 2 | fn / 4 … fn / 2 | 8 |
| 1 | fn / 2 … fn | 16 |
Рисунок 4 - Представление ДВП в частотной области
Вейвлет-фильтрация состоит из следующих действий:
1. прямое дискретное вейвлет преобразование;
2. пороговая обработка коэффициентов, связанная с оценкой средних значений вейвлет коэффициентов на каждом из уровней разложения и последующим пороговым обнулением коэффициентов отражающих слабые импульсные свойства сигнала;
3. обратное вейвлет преобразование.
На рисунке ниже приведем предполагаемую схему организации работы программы.
Рисунок 5 — Организация Вейвлет преобразования.
На рисунке 5 мы видим модуль основной программы, модуль библиотеки gsl, модуль, реализующий параллельное вейвлет преобразование.
Основная программа формирует сигнал, подготавливает материнский вейвлет с помощью модуля GSL и выполняет параллельное Вейвлет-преобразование с помощью модуля cudaDwt.
Рисунок 6 – Ход вейвлет преобразования.Исходя из целей, необходимо разработать модуль, который можно будет использовать без перекомпиляции на платформах одного типа и обладающий переносимостью на уровне компиляции при переносе с платформы на платформу (например, linux/Windows).
Опишем характеристики каждого из модулей:
Модуль CudaDwt — исполняемая статическая/динамическая библиотека, переносимость на уровне исходных кодов между платформами ОС.
Модуль GSL — свободно распространяемая библиотека для научных вычислений, переносимость на уровне исходных кодов между платформами ОС.
Требования:
CPU: intel-совместимый универсальный процессор не ниже i386.
GPU: NVIDIA видеокарта с поддержкой CUDA начиная с 8000-ой серии.
ОС: Windows XP и выше, Linux, Mac OS.
Вейвлет преобразование выполняет свертку сигнала с фильтрами. Что и обуславливает длительную обработку при последовательной реализации.
Рассмотрим подробнее каждый из этапов.
Прямое преобразование
Выполняется сверка целевого сигнала с низкочастотным и высокочастотными фильтрами с одновременным прореживанием, что при последовательном выполнении дает трудоемкость кратную N/2.
Но так как данную операцию на исходной последовательности данных нужно проводить Log2(N) раз, каждый раз уменьшая размерность на 2, то получаем трудоемкость:
Трешолдинг
При пороговой обработке необходимо дважды пройти полученную свертку с низкочастотным и высокочастотными фильтрами и произвести отсечение по уровням. Таким образом, имеем трудоемкость:
Обратное преобразование.
Обратное преобразование производится аналогично прямому и имеет ту же трудоемкость:
При использовании любых последовательных средств вычислений мы сталкиваемся с жесткой зависимостью от объема входных данных, выраженную в трудоемкости:
1.3 Фильтрация на GPGPU
GPGPU – сопроцессор к центральному процессору для массово параллельных вычислений. Теоретически при использовании данного устройства для алгоритмов, допускающих мелкозернистый параллелизм, будет отсутствовать сильная связь с объемом входных данных. То есть при любом объеме мы будем получать результат за ограниченное время.
Перечислим возможности GPGPU как массивно параллельного устройства вычисления:
· множество процессоров виртуально может быть представлено в любой конфигурации с тремя с измерениями;
· каждый процессор имеет уникальный идентификатор;
· процессоры разделены на сильно связные компоненты – потоки в блоках, и слабосвязные – блоки в сетке;
· для оптимизации работы реализована иерархия памяти, технологии группового чтения данных, кэшированного чтения данных;
· одновременно можно запустить на выполнение множество тысяч потоков.
Представим каждую из трех частей в параллельном виде, доступном для запуска на массивно параллельном вычислителе.
Прямое вейвлет преобразование.
При прямом вейвлет преобразовании выполняется свертка и для каждого значения ряда сигнала вычисляется новое значение в соответствии с массивами коэффициентов.
Так как значение свертки для отсчета не зависит от вычисления других, то все вычисления можно распараллелить по номеру отсчета.
Поэтому всё что нужно, это, используя необходимые CUDA расширения языка Си, переписать исходный код последовательной реализации, убрать конструкцию цикла, и сделать вычисление адреса текущего элемента во входной последовательности через уникальные идентификаторы процессора потока.
Что и выполнено ниже:
__global__ void
cudaDWTStepForward(float *a, const size_t n, const size_t nmod, const size_t stride, const float * const h1, const float* const g1, const size_t nc, float*scratch)
{
const size_t idx = blockDim.x * blockIdx.x + tidx; // индекс элемента данных в первой половине = РазмерБлока*номерБлока + индексПотокаВБлоке
const size_t n1 = n - 1;
const size_t nh = n >> 1;
if(idx < nh){
float h = 0, g = 0;
const size_t ni = 2*idx + nmod;
for (size_t k = 0; k < nc; k++) // работа с коэффициентами
{
const size_t jf = n1 & (ni + k);
const float ani = a[stride * jf];
h += h1[k] * ani;
g += g1[k] * ani;
}
scratch[idx] += h; // Преобразованные данные
scratch[idx + nh] += g;
}
Также необходимо произвести расчет числа процессоров под задачу прямого преобразования исходя из формулы:
Число процессоров = числу входных элементов.
Вызов функции на выполнение:
// параметры конфигурации исполнения
const size_t block_size = BL_SZ;
dim3 threads( block_size, 1);
const size_t grid_size = (n/2) / (block_size - nc) + (((n/2)%(block_size - nc))>0); // n/2 потому проходится лишь первая половина
dim3 grid(grid_size, 1);
// запуск ядра на видео-карте
cudaDWTStepForward<<< grid, threads >>>(d_a, n, nmod/*, stride*/, d_h1, d_g1, nc, d_scratch);
Теоретическая трудоемкость данного кода будет приближаться к O(Log2(N)), так как необходимо последовательно выполнить данную функцию для прямого преобразования в соответствии с алгоритмом, описанном выше.
Пороговая обработка
При пороговой обработке мы получаем аналогичную ситуацию и для каждого элемента данных можно запустить свой поток обработки, что продемонстрировано ниже:
__global__ void
cudaDenoiseKernel(float *data, const size_t n,const float threshold)
{
const size_t idx = blockDim.x * blockIdx.x + threadIdx.x; // индекс элемента данных в первой половине = РазмерБлока*номерБлока + индексПотокаВБлоке
if(idx < n){
float d = data[idx];
const float fd = fabsf(d);
float t = fd - threshold;
t = (t + fabsf(t)) / 2.f;
// Signum
if (d != 0.f) d = d / fd * t;
data[idx] = d;
}
}
Для успешного выполнения обработки необходимо жестко контролировать нахождение алгоритма внутри допустимого множества вычисляемых данных. Контроль необходимо осуществлять как при проектировании ядра, так и во время вычисления количества и конфигурации потоков для запуска ядер:
// параметры конфигурации исполнения
const size_t block_size = 16;
dim3 threads( block_size, 1);
const size_t grid_size = len / (block_size) + ((len%(block_size))>0); // n/2 потому проходится лишь первая половина
dim3 grid(grid_size, 1);
// запуск ядра на видео-карте
cudaDenoiseKernel<<< grid, threads >>>(d_a, len, threshold);
Обратное вейвлет преобразование.
Обратное вейвлет преобразование мало отличается от прямого, разве что формулой вычисления.
Также для каждого входного элемента можно создать отдельный поток.
Рассмотренные выше методы распараллеливания довольно тривиальны и очевидны, но, к сожалению, использование только лишь здравого смысла как показывает практика - не дает роста производительности по сравнению с реализацией на универсальном вычислительном устройстве. А дает лишь набор проблем при работе с модулем и неточности результатов.
1.4 Проблемы распараллеливания
При распараллеливании задачи в лоб, без использования дополнительных знаний архитектуры и методов разработки под массивно параллельный вычислители сложно получить эффективную программу.
Приведем проблемы, которые необходимо преодолеть для получения эффективной программы.
· Проблема точности вычислений: GPGPU имеют короткую историю и фактически начали активно развиваться лишь с 2007-го года: в настоящее время поддержка плавающих чисел двойной точности всё еще является предметом обсуждения и встроена не во все имеющиеся GPGPU; также следует отметить потерю производительности при использовании плавающих чисел двойной точности.