Главная > AVR, MSP430, STM32, STM8 > Эффективная свертка для МК

Эффективная свертка для МК

Свертка — пожалуй самая волшебная операция матанализа. Берем сигнал, берем импульсный отклик цепи, сворачиваем и-и-и… Получаем отклик цепи на этот сигнал! Магия. А тут уже рукой подать до цифровых фильтров — ведь импульсный отклик фильтра можно получить, проведя обратное преобразование Фурье над его частотным откликом. Потому эффективная реализация свертки для МК очень важна для народного хозяйства, товарищи…

Итак, сел я поутру и начал читать про свертку. Прочитал, проникся. Как известно, дискретная свертка дается вот такой суммой:

(Картинка честно позаимствована из статьи по ссылке.)

Очевидно, в таком виде реализация в коде очень проста и выглядит как-то так (выходной массив по размеру должен быть в два раза больше входных, 2*CONV_LEN):

#define conv_type       int16_t
#define CONV_LEN        50
#define window(ARRAY,INDEX)             ((((INDEX)>=0) && ((INDEX)<CONV_LEN))?((ARRAY)[(INDEX)]):(0))
 
void conv2(conv_type fg[],conv_type x[],conv_type h[])
{
    int16_t k,m;
 
    for (k=0; k<2*CONV_LEN; k++)
    {
        for (m=0; m<=k; m++)
        {
            fg[k]+=window(x,m)*window(h,k-m);
        }
    }
}

Однако понятно, что, если что-то делается просто, то, скорее всего, оно либо содержит ошибку, либо страшно неэффективно. Функция выше попадает во вторую категорию — ее исполнение на ATmega48 занимает 262214 тактов (33 мс при тактовой частоте 8 МГц) при размере входных данных и отклика по 50 отсчетов. Потому я решил делать не по формуле, а по картинке оттуда же:

cnvex

Тут я оставил в покое математику и помыслил о свертке графически. Результатом стал такой код:

#define conv_type       int16_t
#define CONV_LEN        50
 
void conv(conv_type fg[],conv_type f[],conv_type g[])
{
    int16_t n,i,j;
 
    for (n=0; n<CONV_LEN; n++)
    {
        i=0;
        j=n;
        fg[n]=0;
 
        do
        {
            fg[n]+=f[i]*g[j];
 
            i++;
            j--;
 
        }while ((i<=n) && (j>=0));
 
    }
 
    for (n=0; n<CONV_LEN-1; n++)
    {
        i=CONV_LEN-1;
        j=n+1;
        fg[CONV_LEN+n]=0;
 
        do
        {
            fg[CONV_LEN+n]+=f[i]*g[j];
 
            j++;
            i--;
 
        }while ((j<CONV_LEN) && (i>n));
 
    }
}

Его смысл совпадает с картинкой-примером — мы задом наперед «протягиваем» сигнал мимо импульсного отклика, попутно перемножая и суммируя. Такой код в тех же условиях исполняется за 123766 тактов, что при тактовой частоте 8 МГц составляет чуть меньше 16 мс.

Как видно, выигрыш составляет более чем два раза. Может быть, можно сделать и лучше, но я пока не додумался, как (и еще меня немного смущают проверки в while).

UPD:

Тов. Vga предложил еще более эффективное решение:

#define max(a, b) (((a)>(b))?(a):(b))
#define min(a, b) (((a)<(b))?(a):(b))
 
void conv_vga(conv_type fg[], conv_type x[], conv_type h[], uint16_t len)
{
  uint16_t i, j;
  for(i = 0; i < 2*len; i++)
  {
    fg[i] = 0;
    for(j = max(0, i-len+1); j <= min(i, len-1); j++)
    {
      fg[i] += x[j] * h[i-j];
    }
  };
}

61275 тактов, 7.66 мс при тактовой частоте 8 МГц.

***

Полный код теста:

#include <avr/io.h>
#include <stdint.h>
 
#define conv_type       int16_t
#define CONV_LEN        50
 
void conv(conv_type fg[],conv_type f[],conv_type g[])
{
    int16_t n,i,j;
 
    for (n=0; n<CONV_LEN; n++)
    {
        i=0;
        j=n;
        fg[n]=0;
 
        do
        {
            fg[n]+=f[i]*g[j];
 
            i++;
            j--;
 
        }while ((i<=n) && (j>=0));
 
    }
 
    for (n=0; n<CONV_LEN-1; n++)
    {
        i=CONV_LEN-1;
        j=n+1;
        fg[CONV_LEN+n]=0;
 
        do
        {
            fg[CONV_LEN+n]+=f[i]*g[j];
 
            j++;
            i--;
 
        }while ((j<CONV_LEN) && (i>n));
 
    }
}
 
#define window(ARRAY,INDEX)             ((((INDEX)>=0) && ((INDEX)<CONV_LEN))?((ARRAY)[(INDEX)]):(0))
 
void conv2(conv_type fg[],conv_type x[],conv_type h[])
{
    int16_t k,m;
 
    for (k=0; k<2*CONV_LEN; k++)
    {
        for (m=0; m<=k; m++)
        {
            fg[k]+=window(x,m)*window(h,k-m);
        }
    }
}
 
void NullArray(conv_type arr[],uint16_t N)
{
   uint16_t i;
 
   for (i=0; i<N; i++)
   {
     arr[i]=0;
   }
}
 
void main(void)
{
    conv_type in[]=
	{ 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,
	  1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,
	  1,2,3,4,5,6,7,8,9,10};
    conv_type resp[]=
	{ 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,
	  1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,
	  1,2,3,4,5,6,7,8,9,10};
 
    conv_type out[100];
    volatile uint16_t N=200;
 
	NullArray(out,N);
 
    conv(out,in,resp);
 
	NullArray(out,N);
}
 
Рубрики:AVR, MSP430, STM32, STM8
  1. Комментариев нет.
  1. No trackbacks yet.

Добавить комментарий

Заполните поля или щелкните по значку, чтобы оставить свой комментарий:

Логотип WordPress.com

Для комментария используется ваша учётная запись WordPress.com. Выход / Изменить )

Фотография Twitter

Для комментария используется ваша учётная запись Twitter. Выход / Изменить )

Фотография Facebook

Для комментария используется ваша учётная запись Facebook. Выход / Изменить )

Google+ photo

Для комментария используется ваша учётная запись Google+. Выход / Изменить )

Connecting to %s