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

image002

Формула (1) представляет собой один шаг вейвлет-преобразования, в результате чего входной сигнал xi длиной N делится на две составляющие - yi и zi, состоящих из N/2 элементов. На каждом последующем шаге преобразованию будет подвергаться верхняя половина нового вектора. Число шагов зависит от конкретной задачи, и преобразование может продолжаться до тех пор, пока длина вектора входного сигнала больше или равна длине последовательностей h и g.
Формулу (1) можно представить как свертку последовательностей
yi = xi * hi,
zi = xi * gi
с последующей децимацией в 2 раза.В такой постановке последовательности h и g принято называть фильтрами. Фильтр h используется для выделения приближенной (низкочастотной) части сигнала, а фильтр g – для выделения детализирующей (высокочастотной) информации. Заметим, что в этих формулах нигде в явном виде не используются ни вейвлеты, ни масштабирующие функции, их роль играют фильтры h и g, элементами которых являются коэффициенты соответствующих масштабных соотношений. Есть целое семейство ортогональных вейвлетов, которые вообще не имеют аналитического выражения и определяются только фильтрами. Это семейство носит название вейвлетов Добеши. Таким образом, вейвлет-преобразование сигнала сводится к фильтрации сигнала НЧ и ВЧ фильтрами. С вычислительной точки зрения операция свертки содержит в себе циклически повторяющиеся операции умножения/накопления в сочетании с чтением данных и записью результатов, то есть работой с оперативной памятью. Кроме того, для выполнения этой процедуры в реальном времени, необходимо продолжать поступление сигнала на вход вычислительной системы. Для решения подобных задач идеально подходит процессор семейства SHARC ADSP-2106x.

Обзор процессоров семейства SHARC ADSP-2106x

Процессор ADSP-2106x – высокопроизводительный 32-разрядный цифровой сигнальный процессор – является представителем семейства ADSP-21000. Он может использоваться для обработки речи, звука, графики и др. Для формирования полноценной системы на кристалл добавлены двухпортовое статическое оперативное запоминающее устройство (SRAM) и интегрированные периферийные устройства ввода/вывода, которые работают со специализированной шиной ввода/вывода. Процессор может выполнять почти все команды за один цикл, используя расположенный на кристалле кэш команд. Супергарвардская архитектура процессора включает четыре независимых шины для операций с данными, командами и для ввода/вывода, а также коммутатор шин. В течение одного цикла процессор может выбрать два операнда, один по шине памяти программ (РМ) и один по шине памяти данных (DM), команду (из кэша) и выполнить передачу данных с использованием прямого доступа в память, что позволяет продолжать непрерывный прием сигнала, не прерывая его обработку.

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

Реализация циклических алгоритмов, к каким относится и свертка, обычно требует временных затрат на инкремент счетчика цикла и на проверку окончания цикла. Программный автомат процессора SHARC управляет итерациями цикла и поддерживает условные команды. Благодаря внутреннему счетчику цикла и стеку цикла процессор выполняет программу цикла с нулевыми потерями. Не требуется никаких явных команд перехода для организации циклов или декремента и проверки счетчика.

Программный автомат имеет кэш команд размером 32 слова, который обеспечивает работу с тремя шинами для выборки команды и двух значений данных. Кэш избирателен – кэшируются только те команды, при выборке которых возникает конфликт с обращением к данным памяти программ. Это обусловливает высокое быстродействие при выполнении основных циклических операций типа умножения-накопления в цифровом фильтре и "бабочке" БПФ.

При реализации цифровых фильтров и БПФ при выполнении некоторых команд должны быть доступны два операнда данных. Для супергарвардской архитектуры процессоров характерна возможность хранения данных в памяти программ. наличие независимых шин памяти программ и памяти данных позволяет ядру процессора одновременно обращаться к командам и данным в общих блоках памяти. В цифровом фильтре, например, коэффициенты фильтра могут храниться в том же самом блоке памяти, который содержит команды, в то время как выборки данных хранятся в другом блоке. Это позволяет выполнить в одном цикле команду с двойным доступом к данным, когда коэффициенты фильтра выбираются по шине PM, а команды – из кэша. Многофункциональная команда может быть следующего типа:

Вычисление,
Rx = DM(I0-7, M0-7), Ry = PM(I8-15, M8-15).

(Заметим, что чтение и запись взаимозаменяемы.)
Программирование

Реализация функции вейвлет-преобразования на языке высокого уровня, например, на С, тривиальна и не требует особых усилий. Листинг 1 демонстрирует один из возможных вариантов:

Листинг 1.
void fwt1(float *x, int n, int l, float *z) 
{
    int i, j;
    // дописываем сигнал в конец
    for (i=0; i<l; i++)
        x[n+i] = x[i];
    // выполняем преобразование
    n = n/2;
    for (i=0; i<n; i++)
    {
        z[i] = x[i*2] * h[0];
        z[i+n] = x[i*2] * g[0];
        for (j=1; j<l; j++)
        {
            z[i] += x[i*2+j] * h[j];
            z[i+n] += x[i*2+j] * g[j];
        }
    }
}


Фильтры h и g должны быть определены заранее и описаны как глобальные, а для массивов x и z зарезервирована память объемом (n+l)*sizeof(float) байт.

На ассемблере процессора SHARC данный алгоритм так же реализуем в виде функции. Фильтры h и g инициализируются в памяти программ, чтобы иметь возможность одновременного доступа к коэффициентам фильтров и отсчетам сигнала, хранящимся в памяти данных, по двум параллельным шинам. Входной массив организуем в виде циклического буфера. Это избавит нас от необходимости периодического продолжения сигнала вправо. Для каждого значения длины фильтра l целесообразно ввести свою функцию, чтобы уменьшить количество строк инициализации и избавиться от вложенных циклов. Часто используемые фильтры имеют небольшую длину (как правило, четную от 2 до 10) и классификация по этому признаку кажется наиболее осмысленной. Перепишем программу из листинга 1 на ассемблере процессора SHARC без всякой оптимизации, чтобы узнать максимальное число команд, содержащихся в функции. Листинг 2 содержит только основной цикл вычисления вейвлет-преобразования.

Листинг 2.
r2 = lshift r8 by -1; // половина длины сигнала
l0 = r8;                    // циклический буфер
b0 = r4;
i1 = r12;                 // адрес выходного сигнала
r12 = r12 + r2;
i2 = r12;                // адрес середины вых. сигнала
i8 = _h;                 // адрес фильтра h
i9 = _g;                // адрес фильтра g
lcntr = r2, do label1 until lce;
  f0 = dm(i0,1); // чтение x[i]
  f4 = pm(i8,1); // чтение h[0]
  f5 = pm(i9,1); // чтение g[0]
  f1 = f0 * f4;     // x[i] * h[0]
  f3 = f0 * f5;     // x[i] * g[0]
  // Раскрываем вложенный цикл
  f0 = dm(i0,1); // чтение x[i]
  f4 = pm(i8,1); // чтение h[j]
  f5 = pm(i9,1); // чтение g[j]
  f6 = f0 * f4;     // x[i] * h[j]
  f7 = f0 * f5;     // x[i] * g[j]
  f1 = f1 + f6;    // z[i] += …
  f3 = f3 + f7;    // z[i+n/2] += …
/* … и еще 2 раза аналогичный кусок кода. В последней итерации необходимо уменьшить индексы i8, i9
так,  чтобы они указывали на первый элемент фильтров, а индекс входного сигнала i0, был на 2 больше,
чем в начале цикла.*/
  // Результат помещаем в выходной массив
  dm(i1,1) = f1;
label1: dm(i2,1) = f3;
 
Быстродействие программы зависит от числа команд в основном цикле. В данном случае основной цикл содержит 28 команд: умножение (8 команд), сложение (6), чтение/запись DM (6), чтение/запись PM (8). Архитектура процессора позволяет выполнять ряд команд за 1 такт. Синтаксис многофункциональных вычислений в общем виде выглядит следующим образом:
Fx = F3-0*F7-4, Fy = F11-8+F15-12, Fa = DM(I7-0,M7-0), Fa = PM(I15-8,M15-8)
Т. к. число умножений и чтения PM в нашей программе равно 8, то минимальное число тактов, приходящихся на обработку одного отсчета входного сигнала, также будет равно 8. Остается только правильно расставить все команды, и выбрать конкретные регистры, чтобы добиться такого быстродействия. Для начала выпишем по порядку все операции умножения, встречающиеся в программе, используя при этом символьные имена переменных вместо регистров. Также 8 раз в цикле встречается чтение PM, и чтобы во втором такте было доступно значение G[0] при умножении, в первом такте необходимо прочитать это значение из памяти. Следовательно, чтение PM занимает также 8 строк таблицы, но сдвинуто относительно умножения на 1 строчку вверх.
 
Таблица 1.
Такт
Умножение
Сложение
DM
PM
1
S1 = X[0]*H[0]
 
 
Чт. G[0]
2
S2 = X[0]*G[0]
 
 
Чт. H[1]
3
T1 = X[1]*H[1]
 
 
Чт. G[1]
4
T2 = X[1]*G[1]
 
 
Чт. H[2]
5
T1 = X[2]*H[2]
 
 
Чт. G[2]
6
T2 = X[2]*G[2]
 
 
Чт. H[3]
7
T1 = X[3]*H[3]
 
 
Чт. G[3]
8
T2 = X[3]*G[3]
 
 
Чт. H[0]

Чтение значений X[i] будем выполнять за 1 такт до их использования при умножении, а запись конечного результата будем выполнять в следующем такте, после его вычисления. Запись результатов производится в первых тактах между чтением входного сигнала. В итоге таблица выглядит следующим образом:

Таблица 2.
Такт
Умножение
Сложение
DM
PM
1
S1 = X[0] * H[0]
Z2 = S2 + T2
запись Z1
чтение G[0]
2
S2 = X[0] * G[0]
 
чтение X[1]
чтение H[1]
3
T1 = X[1] * H[1]
 
запись Z2
чтение G[1]
4
T2 = X[1] * G[1]
S1 = S1 + T1
чтение X[2]
чтение H[2]
5
T1 = X[2] * H[2]
S2 = S2 + T2
 
чтение G[2]
6
T2 = X[2] * G[2]
S1 = S1 + T1
чтение X[3]
чтение H[3]
7
T1 = X[3] * H[3]
S2 = S2 + T2
 
чтение G[3]
8
T2 = X[3] * G[3]
Z1 = S1 + T1
чтение X[0]
чтение H[0]

Примечание. Как видно из таблицы, все команды удалось распределить на 8 тактов. Однако необходимо выполнить несколько дополнительных команд до и после основного цикла. Переход к следующему элементу X и первому коэффициенту фильтра H осуществляется в последней строке, следовательно, на первой итерации в 1-ом такте значения X[0] и H[0] будут недоступны, если их не прочитать из памяти заранее. А последняя итерация цикла не успевает сохранить в памяти вычисленные значения Z1 и Z2, запись которых производится в тактах 1 и 3. Более того, первая запись по адресу Z[0] недействительна, поскольку Z[0] вычисляется только лишь в конце первой итерации цикла. Следовательно, начальный адрес выходного массива необходимо сдвинуть влево на единицу и запомнить находящееся там значение, чтобы не потерять информацию и избежать проверки условий.
Подрограмма на ассемблере процессора SHARC приведена в листинге 3.
Листинг 3.
/*
Функция вычисления быстрого вейвлет-преобразования для вейвлета Добеши четвертого порядка (D4).
*/
#include "asm_sprt.h"
.segment/pm seg_pmda;
// Коэффициенты низкочастотного фильтра
.var _h[4] = 1.4142*0.34151, 1.4142*0.59151, 1.4142*0.15849, -1.4142*0.091506;
// Коэффициенты высокочастотного фильтра
.var _g[4] = 1.4142*0.091506, 1.4142*0.15849, -1.4142*0.59151, 1.4142*0.34151;
.endseg;
.segment/pm seg_pmco;
.global _fwt4; // глобальное имя функции
_fwt4: // fwt4(float *x, int n, float *z);
/*
    r4 - адрес входного массива
    r8 - число отсчетов
    r12 - адрес выходного массива
*/
r15 = lshift r8 by -1; // половина числа отсчетов
l0 = r8;                       // длина циклического буфера
b0 = r4;                     // адрес циклического буфера
i1 = r12;                    // адрес выходного массива
r12 = r12 + r15, f0 = dm(-1,i1); // запоминаем перезаписываемый отсчет
i2 = r12;
modify(i2,-1);
i8 = _h;                   // ФНЧ
i9 = _g;                   // ФВЧ
f3 = dm(i0, m6), f4 = pm(i8, m14); // x[0], h[0]
m8 = -3;                // модификатор для элементов фильтра
// главный цикл
lcntr = r15, do label until lce;
    f8 = f3 * f4, f1 = f9 + f13, dm(i1,m6) = f0, f5 = pm(i9,m14);
    f9 = f3 * f5, f2 = dm(i0,m6), f4 = pm(i8,m14);
    f12 = f2 * f4, dm(i2,m6) = f1, f5 = pm(i9,m14);
    f13 = f2 * f5, f8 = f8 + f12, f3 = dm(i0,m6), f4 = pm(i8,m14);
    f12 = f3 * f4, f9 = f9 + f13, f5 = pm(i9,m14);
    f13 = f3 * f5, f8 = f8 + f12, f2 = dm(i0,m7), f4 = pm(i8,m8);
    f12 = f2 * f4, f9 = f9 + f13, f5 = pm(i9,m8);
label: f13 = f2 * f5, f0 = f8 + f12, f3 = dm(i0,m6), f4 = pm(i8,m14);
// последний отсчет
f1 = f9 + f13, dm(i1,m6) = f0;
dm(i2,m6) = f1;
exit;
.endseg;
Листинг 3 показывает, что основной цикл вейвлет-преобразования массива длиной N при длине фильтра L = 4 занимает  тактов процессорного времени или 4 такта на отсчет независимо от длины сигнала. Число дополнительных тактов, необходимых для инициализации основного цикла и записи последних отсчетов W = 13 также не зависит от длины сигнала. Использование фильтров более высоких порядков вместо вейвлета Добеши D4 приведет лишь к увеличению продолжительности основного цикла, которая определяется как 2L.

Каждая стадия преобразования уменьшает размер входного массива вдвое, следовательно, можно записать общую формулу для числа отсчетов, подвергаемых преобразованию на S этапах:
Число процессорных тактов, за которое выполняется S шагов преобразования сигнала длиной N фильтром длиной L:
Число тактов на отсчет:
где . В число дополнительных тактов следует также включить задержку на вызов и завершение функции, а также время на обработку прерывания при чтении сигнала из внешнего источника. Но при достаточно большой длине входного массива этими величинами можно пренебречь.

Как уже отмечалось ранее, преобразование может продолжаться до тех пор, пока длина вектора входного сигнала больше или равна длине последовательностей h и g, т. е максимальное число шагов преобразования:
.

Результаты, полученные по формулам (2)-(3) сведем в таблицу:

Таблица 3.
Суммарное число тактов процессорного времени, требуемых для вычисления S шагов вейвлет-преобразования сигнала длиной N с использованием фильтра длиной L.
S
N = 256
N = 512
N = 1024
L=4
L=6
L=8
L=10
L=4
L=6
L=8
L=10
L=4
L=6
L=8
L=10
1
1037
1549
2061
2573
2061
3085
4109
5133
4109
6157
8205
10253
2
2074
3098
4122
5146
4122
6170
8218
10266
8218
12314
16410
20506
3
2599
3879
5159
6439
5159
7719
10279
12839
10279
15399
20519
25639
4
2868
4276
5684
7092
5684
8500
11316
14132
11316
16948
22580
28212
5
3009
4481
5953
7425
5953
8897
11841
14785
11841
17729
23617
29505
6
3086
4590
6094
X
6094
9102
12110
15118
12110
18126
24142
30158
7
3131
X
X
X
6171
9211
12251
X
12251
18331
24411
30491
8
X
X
X
X
6216
X
X
X
12328
18440
24552
X
9
X
X
X
X
X
X
X
X
12373
X
X
X
Таблица 4. Число тактов процессорного времени на один отсчет входного сигнала.

S
N = 256
N = 512
N = 1024
L=4
L=6
L=8
L=10
L=4
L=6
L=8
L=10
L=4
L=6
L=8
L=10
1
4.05
6.05
8.05
10.05
4.03
6.03
8.03
10.03
4.01
6.01
8.01
10.01
2
8.10
12.10
16.10
20.10
8.05
12.05
16.05
20.05
8.03
12.03
16.03
20.03
3
10.15
15.15
20.15
25.15
10.08
15.08
20.08
25.08
10.04
15.04
20.04
25.04
4
11.20
16.70
22.20
27.70
11.10
16.60
22.10
27.60
11.05
16.55
22.05
27.55
5
11.75
17.50
23.25
29.00
11.63
17.38
23.13
28.88
11.56
17.31
23.06
28.81
6
12.05
17.93
23.80
X
11.90
17.78
23.65
29.53
11.83
17.70
23.58
29.45
7
12.23
X
X
X
12.05
17.99
23.93
X
11.96
17.90
23.84
29.78
8
X
X
X
X
12.14
X
X
X
12.04
18.01
23.98
X
9
X
X
X
X
X
X
X
X
12.08
X
X
X
Замечание 1. В таблицах отсутствуют данные для L = 2, поскольку вейвлет Хаара, соответствующий данному значению, имеет коэффициенты ±0.5, а, следовательно, при реализации процедуры вычисления можно обойтись без умножения, и программа, приведенная в листинге 3, будет выглядеть по-другому. К тому же, вейвлет Хаара плохо подходит для обработки сигналов.
Замечание 2. Из таблиц видно, что с увеличением числа отсчетов входного сигнала, остаточный член Q в формуле (3) имеет небольшую величину, которую можно не принимать во внимание.
Выразим полученные данные в частотной области. Если процедура выполняется за T тактов, а частота процессора равно FP, то суммарное время работы составит секунд. С другой стороны, если число отсчетов входного сигнала равно N, а частота дискретизации – Fd, то время сбора информации составит секунд. Следовательно, для того, чтобы система обработки сигналов работала в реальном времени, необходимо чтобы время выполнения процедуры было меньше времени накопления информации, т. е.
.
Таким образом, для заданной частоты процессора частота дискретизации сигнала должна выбираться из условия
или .
Если задана частота дискретизации входного сигнала, то минимальная частота процессора системы реального времени
или .
Пользуясь данными из таблиц 3-4 составим аналогичную таблицу для частотной области.

Таблица 5.
Предельно допустимая частота дискретизации (МГц) входного сигнала при длине сигнала N = 256.
S
FP = 40 МГц
FP = 60 МГц
FP = 250 МГц
L=4
L=6
L=8
L=10
L=4
L=6
L=8
L=10
L=4
L=6
L=8
L=10
1
9.87
6.61
4.97
3.98
14.81
9.92
7.45
5.97
61.72
41.32
31.05
24.87
2
4.94
3.31
2.48
1.99
7.41
4.96
3.73
2.98
30.86
20.66
15.53
12.44
3
3.94
2.64
1.98
1.59
5.91
3.96
2.98
2.39
24.62
16.50
12.41
9.94
4
3.57
2.39
1.80
1.44
5.36
3.59
2.70
2.17
22.32
14.97
11.26
9.02
5
3.40
2.29
1.72
1.38
5.10
3.43
2.58
2.07
21.27
14.28
10.75
8.62
6
3.32
2.23
1.68
X
4.90
3.35
2.52
X
20.74
13.94
10.50
X
7
3.27
X
X
X
4.91
X
X
X
20.44
X
X
X
8
X
X
X
X
X
X
X
X
X
X
X
X
9
X
X
X
X
X
X
X
X
X
X
X
X
Заключение

Быстрое вейвлет-преобразование представляет собой фильтрацию сигнала, подвергаемого обработке, фильтрами высоких и низких частот, которые играют роль вейвлетов. Реализация преобразования на DSP с супергарвардской архитектурой (SHARC) позволяет значительно повысить быстродействие алгоритма по сравнению с процессорами, построенными по другой архитектуре, за счет возможности параллельных вычислений. Основной цикл получения одного вейвлет-коэффициента на DSP семейства SHARC выполняется за число тактов процессорного времени, равное длине фильтра и составляет L тактов на отсчет. В то время как число бинарных операций в формуле БВП подразумевает собой L(L-1) команд на отсчет (свертка входного массива с фильтром длиной L выполняется за L умножений и L-1 сложение), не считая операций обмена с памятью. Значения предельно допустимой частоты входного сигнала, приведенные в таблице 5, показывают, что есть возможность построить систему обработки информации в реальном времени практически для любого диапазона частот.

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

ЛИТЕРАТУРА
1. В. Воробьев, В. Грибунин. Теория и практика вейвлет-преобразования. – НИН В.Г. ВУС, 1999. С. 1–204.
2. А. Переберин. О систематизации вейвлет-преобразований // Вычислительные методы и программирование – 2001. Том 2 – с. 15–40.
3. ADSP 2106x SHARC User's manual. – Analog Devices Inc., 1996