Лучшие подходы переноса кода MATLAB в фиксированную точку

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

При конвертации проекта из плавающей точки в фиксированную точку инженеры должны определить оптимальные типы данных в фиксированной точке. Эти типы данных должны удовлетворять ограничениям встраиваемой аппаратуры, при этом удовлетворяя системным требованиям по точности вычислений. Fixed-Point Designer помогает разрабатывать алгоритмы в фиксированной точке и конвертировать алгоритмы из плавающей точки в фиксированную точку, автоматически предлагая типы данных и атрибуты арифметики в фиксированной точке. При этом предоставляется возможность сравнения результатов симуляции в фиксированной точке с точностью до бита с эталонными результатами в плавающей точке.

Подготовка кода к переводу в фиксированную точку
Есть три шага, которые следует предпринять для обеспечения плавного процесса конвертации:

1. Отделить основной алгоритм от остального кода.
2. Подготовить код для инструментирования и ускорения.
3. Проверить используемые функции на поддержку фиксированной точки.

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

% ТЕСТОВЫЕ ВХОДЫ
x = randn(100,1);
 
% АЛГОРИТМ
y = zeros(size(x));
y(1) = x(1);
for n=2:length(x)
    y(n) = y(n-1) + x(n);
end
 
% ВЕРИФИКАЦИЯ РЕЗУЛЬТАТОВ
yExpected = cumsum(x);
plot(y-yExpected)
title('Ошибка')

% ТЕСТОВЫЕ ВХОДЫ
x = randn(100,1);
 
% АЛГОРИТМ
y = cumulative_sum(x);
% ВЕРИФИКАЦИЯ РЕЗУЛЬТАТОВ
yExpected = cumsum(x);
plot(y-yExpected)
title('Ошибка')

function y = cumulative_sum(x)
    y = zeros(size(x));
    y(1) = x(1);
    for n=2:length(x)
        y(n) = y(n-1) + x(n);
    end
end

Таблица 1. Код до и после отделения основного алгоритма от тестовой обвязки.

Подготовка алгоритмического кода к инструментированию и ускорению
Инструментирование и ускорение позволяют упростить процесс конвертации. Fixed-Point Designer используется для инструментирования кода и записи минимальных и максимальных значений всех именованных и промежуточных переменных. Этот инструмент может использовать записанные значения, чтобы предложить типы данных для использования в коде в фиксированной точке.

С использованием Fixed-Point Designer можно также ускорять алгоритмы в фиксированной точке путем создания MEX файла, и ускорять симуляции, требуемые для верификации реализации в фиксированной точке относительно оригинальной версии.

Инструментирование и ускорение полагаются на технологию генерации кода, поэтому прежде чем их использовать, требуется подготовить алгоритм к генерации кода – даже если вы не планируете использовать MATLAB Coder или HDL Coder для генерации кода C или кода HDL.

Вначале требуется определить функции или конструкции в вашем коде MATLAB, не поддерживаемые для генерации кода (см. Поддержка языка для получения списка поддерживаемых функций и объектов).

Есть два способа для автоматизации этого шага:

• Добавить директиву

  %#codegen

  в начале файла MATLAB, содержащего ядро вашего алгоритма. Эта директива сообщает анализатору кода, что следует отмечать функции и конструкции, которые не содержатся в подмножестве языка MATLAB, поддерживаемого для генерации кода.
• Использовать инструмент проверки готовности кода для создания отчета, который выявляет обращения к функциям и использование типов данных, не поддерживаемых для генерации кода.

После подготовки алгоритма к генерации кода, можно использовать Fixed-Point Designer для инструментирования и ускорения кода. Используйте

buildInstrumentedMex

для включения инструментирования с целью записи минимальных и максимальных значений всех именованных и промежуточных переменных. Также используйте

showInstrumentationResults

для просмотра отчета о генерации кода с предложенными типами данных в фиксированной точке. Запустите

fiaccel

для трансляции вашего алгоритма MATLAB в MEX файл и ускорения симуляций в фиксированной точке.

Проверка поддержки фиксированной точки функциями, используемыми в алгоритмическом коде

При выявлении функции, не поддерживаемой для генерации кода, у вас есть три возможности:

1. Заменить функцию на эквивалентную функцию в фиксированной точке.
2. Написать собственную эквивалентную функцию.
3. Изолировать неподдерживаемую функцию используя приведение типа к double на входе функции, и обратное приведение типа к фиксированной точке на выходе.

Затем можно продолжить конвертацию кода в фиксированную точку, и вернуться к неподдерживаемой функции, когда у вас будет подходящая замена (Таблица 2).

y = 1/exp(x);

y = 1/exp(double(x));

Таблица 2. Код до и после изоляции операции в плавающей точке при помощи приведения типа (обратное приведение типа к фиксированной точке на выходе не показано).

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

Например, рассмотрим следующий код:

y = y + x(n)

Это выражение перезаписывает y значением

 y + x(n)

Когда вы используете типы данных в фиксированной точке в коде (для y и x), то тип данных y может поменяться после перезаписи, потенциально приводя к росту разрядности.

Для сохранения типа данных y используйте синтаксис

(:) =

(Таблица 3). Такой синтаксис, известный как индексное присваивание, заставляет MATLAB сохранять существующий тип данных и размер массива перезаписываемой переменной. Выражение

y(:) = y + x(n)

приведет выражение справа к оригинальному типу данных y и предотвратит рост разрядности.

y = 0;
for n=1:length(x)
  y = y + x(n);
end

y = 0;
for n=1:length(x)
  y(:) = y + x(n);
end

Таблица 3. Код до и после использования индексного присваивания для предотвращения роста разрядности.

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

Для применения этой лучшей практики надо сделать следующее:

1. Использовать
   

   cast(x,'like',y)

   или
   

   zeros(m,n,'like',y)​

   для приведения типа переменной к желаемому типы данных при первом определении переменной.
2. Создать таблицу определений типов, начиная с оригинальных типов данных, используемых в коде – обычно, в плавающей точке двойной точности – типа данных по умолчанию в MATLAB (Таблица 4а).
3. До конвертации в фиксированную точку, добавить тип данных single в таблицу типов для поиска несоответствий и других проблем (Таблица 4б).
4. Осуществить верификацию привязки, запустив код, привязанный к каждой таблице, с различными типами данных и сравнив результаты.

% Алгоритм
n = 128;
y = zeros(size(x));

% Алгоритм
T = mytypes('double');
n = cast(128,'like',T.n);
y = zeros(size(x),'like',T.y);
 
% Таблица типов
function T = mytypes(dt)
  switch(dt)
    case 'double'
      T.n = double([]);
      T.y = double([]);
   end
end

Таблица 4a. Код до и после создания таблицы типов данных для отделения алгоритмического кода от определений типов.

% Таблица типов
function T = mytypes(dt)
  switch(dt)
    case 'double'
      T.n = double([]);
      T.y = double([]);
   end
end

% Таблица типов
function T = mytypes(dt)
  switch(dt)
    case 'double'
      T.n = double([]);
      T.y = double([]);
      
    case 'single'
      T.n = single([]);
      T.y = single([]);
  end
end

Таблица 4б. Код до и после добавления типа данных единичной точности к таблице типов.

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

Чтобы получить набор предложенных типов данных для вашего кода, используйте команды Fixed-Point Designer

buildInstrumentedMex

и

showInstrumentationResults

(Таблица 5). Вам потребуется набор тестовых векторов, который задействует полный диапазон типов – поскольку типы, предложенные Fixed-Point Designer, хороши настолько, насколько хороши тестовые воздействия. Продолжительная симуляция с широким набором ожидаемых входных данных приведет к лучшим предложенным типам данных. Выберите первоначальный набор данных в фиксированной точке из предложенных в отчете по генерации кода (Рисунок 1).


Рисунок 1. Отчет по генерации кода, созданный showInstrumentationResults с предложенными типами данных для переменных в алгоритме фильтрации.

Затем вы можете подстроить предложенные типы при необходимости (Таблицы 5 и 6).

function [y,z] = myfilter(b,x,z)
  y = zeros(size(x));
  for n=1:length(x)
    z(:) = [x(n); z(1:end-1)];
    y(n) = b * z;
  end
end

% Тестовые входы
b = fir1(11,0.25);
t = linspace(0,10*pi,256)';
x = sin((pi/16)*t.^2);  % Линейный импульсный сигнал
z = zeros(size(b'));
 
% Построение
buildInstrumentedMex myfilter ...
  -args {b,x,z} -histogram
 
% Запуск
[y,z] = myfilter_mex(b,x,z);
 
% Отображение
showInstrumentationResults myfilter_mex ...
  -defaultDT numerictype(1,16) -proposeFL

Таблица 5. Алгоритм фильтрации и тестовый скрипт для инструментирования и выполнения кода и отображения предложенных типов данных в фиксированной точке для переменных.

function [y,z] = myfilter(b,x,z,T)
  y = zeros(size(x),'like',T.y);
  for n=1:length(x)
    z(:) = [x(n); z(1:end-1)];
    y(n) = b * z;
  end
end

% Тестовые входы
b = fir1(11,0.25);
t = linspace(0,10*pi,256)';
x = sin((pi/16)*t.^2);  % Линейный импульсный сигнал
 
% Приведение типов входов
T = mytypes('fixed16');
b = cast(b,'like',T.b);
x = cast(x,'like',T.x);
z = zeros(size(b'),'like',T.x);
 
% Запуск
[y,z] = myfilter(b,x,z,T);

function T = mytypes(dt)
  switch dt
    case 'double'
      T.b = double([]);
      T.x = double([]);
      T.y = double([]);
 
     case 'fixed16'
      T.b = fi([],true,16,15);
      T.x = fi([],true,16,15);
      T.y = fi([],true,16,14);
  end
end

Таблица 6. Тестовый скрипт и алгоритм фильтрации из Таблицы 4 с типами данных в фиксированной точке.

Запустите ваш алгоритм с новыми типами данных в фиксированной точке и сравните выход с результатами эталонного алгоритма в плавающей точке.

Оптимизация типов данных
Выбрали ли вы ваши собственные типы в фиксированной точке или использовали предложенные Fixed-Point Designer, всегда ищите возможности по оптимизации размеров слова, размеров дробной части, знаковости и, возможно, даже режимов арифметики (fimath). Это можно сделать с использованием разрядных чисел двойной точности (scaled doubles), рассматривая гистограмму значений переменной или путем тестирования различных типов данных из вашей таблицы типов.

Использование Scaled Doubles для выявления потенциальных переполнений
Scaled doubles являются гибридами чисел в плавающей точке и фиксированной точке. Fixed-Point Designer хранит scaled doubles в виде чисел двойной точности, но сохраняет информацию о разрядности, знаке и длине слова. Для использования scaled doubles требуется установить свойство перезаписи типа данных (data type override, DTO) (Таблица 7).

numerictype

>> T.a = fi([], 1, 16, 13, 'DataType', 'ScaledDouble');

>> a = cast(pi, 'like', T.a)

a =
    3.1416

  DataTypeMode: Scaled double: binary point scaling
    Signedness: Signed
    WordLength: 16
FractionLength: 13

fipref

>> fipref('DataTypeOverride', 'ScaledDoubles');

>> T.a = fi([], 1, 16, 13);

>> a = cast(pi, 'like', T.a)

a =
    3.1416

  DataTypeMode: Scaled double: binary point scaling
    Signedness: Signed
    WordLength: 16
FractionLength: 13

Таблица 7. Методы установки свойства перезаписи типов данных – локально и глобально.

Не забудьте сбросить глобальный DTO, если в нем больше нет необходимости, при помощи команды

reset(fipref)

Используйте buildInstrumentedMex для запуска вашего кода и showInstrumentationResults для просмотра результатов. В отчете по генерации кода значения, которые бы переполнились, подсвечены красным цветом (Рисунок 2).


Рисунок 2. Отчет по генерации кода, показывающий переполнения при использовании типа Scaled Doubles (слева) и иконка гистограммы (справа).

Проверка распределения значений переменной
Вы можете использовать гистограмму для выявления типов данных со значениями, которые находятся в допустимом диапазоне, вне диапазона или ниже разрешения (точности). Щелкнув по иконке гистограммы, можно запустить NumericTypeScope и увидеть распределение значений, наблюдаемых во время симуляции для выбранной переменной (Рисунок 3).


Рисунок 3. Гистограмма, показывающая распределение значений переменной, которая переполнилась (“Outside range”), показанное красным цветом.

Тестирование различных типов данных из таблицы типов
Вы можете добавлять собственные вариации типов данных в фиксированной точке в таблицу типов (Таблица 8).

function [y,z] = myfilter(b,x,z,T)
  y = zeros(size(x),'like',T.y);
  for n=1:length(x)
    z(:) = [x(n); z(1:end-1)];
    y(n) = b * z;
  end
end

function mytest
  % Тестовые входы
  b = fir1(11,0.25);
  t = linspace(0,10*pi,256)';
  x = sin((pi/16)*t.^2);  % Линейный импульсный сигнал
  
  % Запуск
  y0  = entrypoint('double',b,x);
  y8  = entrypoint('fixed8',b,x);
  y16 = entrypoint('fixed16',b,x);
  
  % Графики
  subplot(3,1,1);plot(t,x,'c',t,y0,'k');
  legend('Вход','Эталонный выход')
  title('Эталонный запуск')
         
  subplot(3,2,3);plot(t,y8,'k');
  title('8-бит фикс. точка: Выход')
  subplot(3,2,4);plot(t,y0-double(y8),'r');
  title('8-бит фикс. точка: Ошибка')
 
  subplot(3,2,5);plot(t,y16,'k');
  title('16-бит фикс. точка: Выход')
  xlabel('Время (с)')
  subplot(3,2,6);plot(t,y0-double(y16),'r');
  title('16-бит фикс. точка: Ошибка')
  xlabel('Время (с)')
end
 
function [y,z] = entrypoint(dt,b,x)
  T = mytypes(dt);
  b = cast(b,'like',T.b);
  x = cast(x,'like',T.x);
  z = zeros(size(b'),'like',T.x);
  [y,z] = myfilter(b,x,z,T);
end

function T = mytypes(dt)
  switch dt
    case 'double'
      T.b = double([]);
      T.x = double([]);
      T.y = double([]);
 
    case 'fixed8'
      T.b = fi([],true,8,7);
      T.x = fi([],true,8,7);
      T.y = fi([],true,8,6);
 
    case 'fixed16'
      T.b = fi([],true,16,15);
      T.x = fi([],true,16,15);
      T.y = fi([],true,16,14);
  end
end

Таблица 8. Тестовый скрипт для изучения эффектов использования различных типов данных в фиксированной точке из таблицы типов для функции фильтрации.

Сравнение результатов разных итераций для верификации точности алгоритма после каждого изменения (Рисунок 4).


Рисунок 4. Графики результатов тестового скрипта из Таблицы 8, показывающие выход и ошибку после конвертации в 8-битные и 16-битные типы данных в фиксированной точке.

Оптимизация алгоритма

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

• Использовать свойства fimath для улучшения эффективности сгенерированного кода
• Заменять встроенные функции более эффективными реализациями в фиксированной точке
• Реализовывать операции деления иными методами

Использование свойств fimath для улучшения эффективности сгенерированного кода

При использовании настроек fimath по умолчанию, может генерироваться дополнительный код для реализации насыщения при переполнении, округления и арифметики с полной точностью (Таблица 9а).

Компилируемый код:

function y = adder(a,b)
  y = a + b;
end

С типами, заданными с настройками fimath по умолчанию:

T.a = fi([],1,16,0);
T.b = fi([],1,16,0);
 
a = cast(0,'like',T.a);
b = cast(0,'like',T.b);

int adder(short a, short b)
{
  int y;
  int i0;
  int i1;
  int i2;
  int i3;
  i0 = a;
  i1 = b;
  if ((i0 & 65536) != 0) {
    i2 = i0 | -65536;
  } else {
    i2 = i0 & 65535;
  }
 
  if ((i1 & 65536) != 0) {
    i3 = i1 | -65536;
  } else {
    i3 = i1 & 65535;
  }
 
  i0 = i2 + i3;
  if ((i0 & 65536) != 0) {
    y = i0 | -65536;
  } else {
    y = i0 & 65535;
  }
 
  return y;
}

Таблица 9а. Оригинальный код MATLAB и код C, сгенерированный с настройками fimath по умолчанию.

Чтобы сгенерированный код был более эффективным, требуется выбрать такие настройки арифметики с фиксированной точкой, которые подходят типам вашего процессора. Используйте свойства fimath для описания арифметики, способов округления и действий при переполнении, чтобы задать правила осуществления арифметических операций с вашими объектами fi (Таблица 9б).

Компилируемый код:

function y = adder(a,b)
  y = a + b;
end

С типами, заданными с настройками fimath по умолчанию:

T.a = fi([],1,16,0); T.b = fi([],1,16,0); a = cast(0,'like',T.a); b = cast(0,'like',T.b);

int adder(short a, short b)
{
  return a + b;
}

Таблица 9б. Оригинальный код MATLAB и код C, сгенерированный с настройками fimath, подходящими к типам процессора.

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

Реализация операций деления другими способами
Операции деления часто не поддерживаются полностью со стороны аппаратуры, и могут привести к медленным вычислениями. Когда в вашем алгоритме требуется операция деления, рассмотрите возможность замены её на более быструю альтернативу. Если знаменатель является степенью двойки, используйте битовый сдвиг; например, используйте bitsra(x,3) вместо x/8. Если знаменатель является константой, умножьте на обратную величину; например, используйте x*0.2 вместо x/5.

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