DSP на OpenCL -- строительные блоки и блочные диаграммы

Но по сути, мы тут не про язык разговариваем. У меня копится материал для изложения, который я выкладываю для обсуждения. Первая часть изложения - про объекты в протоколе и в описании модели. А вторая часть продолжение темы выкручивания алгоритмов DSP для примения в составе блочной диаграммы.

Блоки в языке OpenCL C - это такие чудные вложенные функции со шляпкой. В GCC поддерживается расширение языка Си, т.н. "nested functions", похожие сущности -- гнездящиеся функции. В языке OpenCL C 2.0+ можно на стороне устройства описывать из блоков асинхронный процесс обработки.

Достаточно крупные операции в DSP надо делать inline функциями, вернее давать возможность компилятору оптимизировать(векторизовать) большие блоки кода. Один из основных методов дающих существенный прирост производительности - векторизация циклов. Это значит что одиночные циклы, для которых число итераций определяется снаружи следует писать в виде inline функций. Мне нравится набор функций, который предлагает ARM CMSIS DSP. Но есть оговорки и вопросы оптимизации под аппаратуру. У нас есть задача переносить коэффициенты и модели из настольного компьютера в контроллер. В компьютере мы используем векторизацию float, в устройстве управления и сбора данных, сервере автоматизации используем ARM Neon, планируем использовать OpenCL. Нам нужен объект и формат данных для переноса "тензоров" и организации "потока тензоров" средствами протокола обмена, между платформами, нам нужен единый API (оптимизированный под конкретную архитектуру) для представления блочных диаграмм. Концептуально, мы говорим про внешний интерфейс для представления функционального блока FunctionalBlockObject, объекты-контейнеры данных TensorObject. А также мы будем рассматривать некоторый набор имен программных блоков, представленных в библиотеке DSP, чтобы ссылаться на их использование в блочной диаграмме.
Говоря про формат данных для внешнего представления блочных диаграмм, мы на сегодня будем ориентироваться на форматы принятые в BACnet: правила кодирования ASN.1, формат Json и CSML. Эти форматы мы выбираем прежде всего, чтобы не плодить сущности, поддержка этих форматов уже присутствует в составе нашего ПО.

Тензоры

Тензор - это массив данных определенного типа: float, double, int, char. Мы говорим про внешнее представление данных. Тнезор обладает свойством выравнивания данных - это его внутреннее свойство, определяется автоматически при создании объекта. Тензор имеет порядок укладки данных при записи и чтении, формат записи. Порядок записи включает: число каналов (feature), пакет данных (batch), пространственную размерность (spatial). Например, изображение с камеры имеет размерность 2d (spatial), организовано как поток кадров (batch), может иметь несколько каналов (feature). Звук имеет размерность =1 (spatial), организовано в поток по времени (batch), и разделено на каналы (features). Тензор - это единица обрабатываемой информации, с точки зрения функционального блока DSP или NN.
Мы будем исходить из того, что передать тензор целиком в одном пакете может оказаться не возможно. Поэтому обращение к тензору может быть как к массиву, ReadRange - запрос на чтение фрагмента массива. Мы предполагаем некоторую процедуру заполнения входных данных, объектов связанных между собой блочной диаграммой, затем запуск преобразования. Результатом преобразования может являться набор тензоров. Результат преобразования может распространяться по подписке CoV. Преобразование запускается по готовности входных данных, т.е. когда все входные данные(тензоры) записаны.
Тензоры могут быть вида TensorOutputObject, TensorInputObject. TensorOutputObject - обладает командным свойством на PresentValue. Когда записывается содержимое тензора, происходит запуск зависимостей по готовности.

Функциональные блоки

Функциональный блок содержит ссылку на исполняемый код(kernel) - DeviceObjectReference, список аргументов, список фиксированный идентификаторов объектов ARRAY_OF(DeviceObjectReference) - зависимости. В процессе исполнения, блок оформляет подписку на соответствующие объекты списка зависимостей. Функциональный блок имеет статус(состояние,контекст) - это некоторая структура данных, в которой могут хранится промежуточные величины. Контекст сохраняется и может быть перезаписан в режиме out-of-service.
Для блока мы определяем фазы: Инициализация, обновление контекста, обработка, пауза/завершение/сериализация контекста. Блочная диаграмма это не сам конвейер, это пакетное задание, которое может выполняться на конвейере с использованием данных. Конвейер это очередь исполнения по готовности.
Функциональный блок может обладать отладочной информацией (опция), для отладки используются штампы времени: queued, submit, start, end.

Программы

Программа - это текст или бинарник. Со стороны хост-приложения возможно определить из бинарника список блоков и список агрументов(входных/выходных портов). Каждый порт описывается 'форматом записи', однако формат записи из программы не определить.

Строительные блоки

Список функций API

dot(a,b)

операция свертки вектора sum_{0}^{N-1}(a[i]·b[i])

conv(a,b,n)

операция частичной свертки, выполняется с разворотом вектора sum_{0}^{N-1}(a[i]·b[n-i-1]). В это определение попадает корреляционная функция

float cascade_biquad_df2T(z,x)

вычисление значения каскада y = z+b0·x

float cascade_biquad_df2T_update(S,x)

обновление регистров каскада (S)

float fir(S,x)

FIR фильтр. Интерполяция.

pid(S,x)

PID регулятор.

lms_update(S,x)

LMS подстройка параметров адаптивного фильтра

IIR фильтры обрабатывают batch составляющую, время.
Для задач мультирейт используются фильтры LPF_decimate (по времени).
Аудио процессор может обрабатывать полосы и сшивать их на выходе. Для этого используются фильтры Allpass. В каких-то случаях могут использоваться фильтры iir_lattice.
Пид регулятор - это фильтр типа biquad. LMS фильтр - подстройка.

PID регулятор

 

Передаточная функция:
H(s) = ( K_P + K_I·1/s + K_D·s ) Замена: 1/s = T/2·(1+z^-1)/(1-z^-1), s = 1/T·(1-z^-1)

Передаточная функция в дискретном пространстве:
H(z) = KP + KI·T/2·(1+z-1)/(1-z-1) + KD/T·(1-z-1)
или 
H(z) = KP + KI·T·1/(1-z-1) + KD/T·(1-z-1)
Ki := KI·T
Kd := KD/T
Может быть представлена биквадратной секцией

void pid_init(S, x){
  S->b0 =  Kp + Ki + Kd;

  S->b1 = -Kp -2Kd
  S->b2 =  Kd;

  S->z1 = 0;
  S->z2 = 0;
}
// первый проход
float pid(S, x){
  return b0*x+z1;
}
// второй проход, обновление регистров
void pid_update(S, x){
  z1 =(b1-b0)*x + (z2 - z1);
  z2 = b2*x;
}

При разработке PID регулятора мы также используем два прохода: первый проход встроен в поток обработки, второй - поток обновления регистров и состояния фильтров. Пока что у нас получается что абсолютно все фильтры представлены выражением: y = g*x+Z, включая и PID регулятор. Далее мы хотим добавить свойство автоматической подстройки каэффициентов FIR фильтров и IIR фильтров. На втром цикле обновления регистров, мы добавляем возможность подстройки коэффциентов по методу LMS (). Есть известное преобразование FIR фильтра в LMS - получается адаптивный фильтр. Мы создаем сеть, в которой каждый функциональный блок должен обладать адаптивностью. Точнее можно сказать, что в процессе работы возникают периоды обучения - подстройки коэффициентов модели. Таким образом, фаза UPDATE может сопровождаться "обучением", если известна величина ошибки выходного сигнала. Для каждого слоя сети, мы можем определить этот признак "обучения".

LMS фильтры

Фильтры адаптивные (LMS) мы будем изготавливать из фильтров типа IIR и FIR. Метод LMS можно применять для обучения (подстройки весовых коэффициентов) нейронной сети. Мы вводим единую методику обучения/распознавания для нейронных сетей (NN) и для структур цифровой обработки сигналов (DSP). Нам известен способ изготовления адаптивного цифрового фильтра из FIR фильтра эту методику мы хотим расширить на остальные задачи.

Для некоторого преобразования c передаточной функцией представленной полиномом H(z) = h₀+h₁z^-1+h₂z^-2+...+h_Nz^-N , можно подстраивать коэффициенты, если известна ошибка выходного значения. Ошибка может быть пересчитана в поправку каждого конкретного коэффициента.
Ошибку на выходе определим, как квадрат отклонения.
E = (d[n] - y[n])² = (d[n] - sum(b[i]z^-i))²
изменение коэффицента для фильтра можно определить как градиент. Используем метод "градиентного спуска"
w[n+1]=w[n]-μ·∂E/∂w
Здесь μ - это величина шага подстройки. Для стабильного проявления методики, бывает полезно ограничить ошибку в каких то пределах, используя функцию clamp(). Берем частную производную, не сильно увлекаясь частностями. w_i[n+1] = w_i[n] - 2μ·e[n]·∂e[n]/∂w_i
∂E/∂w_i = -2(e[n])(x[n-i])

b_i[n+1]=b_i[n]+2μ (e[n])(x[n-i]).
Это все! Каждый коэффицент подстраиваем исходя из ошибки на выходе. IIR фильтр - это комбинация двух фильтров FIR. Для подстройки коэффициентов A(z) в IIR-фильтре с передаточной функцией H(z)=B(z)/A(z) используем аналогичный вывод:
a_i[n+1]=a_i[n]+2μ (e[n])(y[n-i]).

{Тут в изложении присутсвует некоторое искусство, с которым надо поразбираться. Работать будет в любом случае, однако математически правильный вариант брать частную производную не от уравнения в конечных разностях. Варианты, которые можно получить исследуя конечные разности и интегралы от дельта-функций - это линии задержки разного вида и фильтры на поток ошибок. К сожалению, быстро ответ не нашел. Нашел, что этим когда-то занимался I.-D. Landau, }
{1976 Feintuch’s algorithm}
{A parallel-type adaptive algorithm is proposed which utilizes the fast least-squares method. Its computational complexity is much less than that of I. Landau's (1978) method, which is based on hyperstability theory Hyperstability requires much computation because it involves matrix operations. The proposed method has nothing to do with hyperstability. It is also shown theoretically and using computer simulation, that the convergence performance of the proposed method is better than that of the series-parallel-type fast least-squares method.} Линия задержки будет разная для разных форм представления фильтров. Глядя на презентацию, мне почему-то кажется, что тут может быть ошибка в интерпретации, что-то не то с частной производной. Должна быть простая форма. Мне видится простая форма при вторичном использовании регистров z. Предположительно выглядит так, доказать надо на практике:

void lms_dt2t_update(S, x, ue){
  b0+= ue*x;
  b += ue*x;// вектор
  a += ue*y;// вектор
}
void lms_update(S, x, ue){
  b0+= ue*x;
  b += ue*z;// вектор
  a += ue*z;// вектор
}

Если мы применяем LMS к фильтру типа PID должны быть веса разные для коэффициентов. Коэффциенты A не подстравиваются.