Самый быстрый метод расчета свертки

Кто-нибудь знает о самом быстром способе вычисления свертки? К сожалению, матрица, с которой я имею дело, очень большая (500x500x200), и если я использую convn в MATLAB это занимает много времени (мне приходится повторять этот расчет во вложенном цикле). Итак, я использовал свертку с FFT и теперь она быстрее. Но я все еще ищу более быстрый метод. Любая идея?

8

Решение

Если ваше ядро ​​отделимо, наибольшее увеличение скорости будет реализовано при выполнении нескольких последовательных одномерных сверток.

Стив Эддинс из MathWorks описывает, как воспользоваться ассоциативностью свертки для ускорения свертки, когда ядро ​​отделимо в контексте MATLAB на его блог. Для P-by-Q ядро, вычислительное преимущество выполнения двух отдельных и последовательных сверток против двумерной свертки PQ/(P+Q), что соответствует 4,5x для ядра 9×9 и ~ 11x для ядра 15×15. РЕДАКТИРОВАТЬ: Интересная невольная демонстрация этой разницы была дана в этот вопрос&.

Чтобы выяснить, является ли ядро ​​отделимым (то есть внешним произведением двух векторов), блог продолжает описывать как проверить, может ли ваше ядро ​​отделиться от SVD и как получить 1D ядро. Их пример для двумерного ядра. Для решения для N-мерной отделимой свертки, проверьте это FEX представление.


Еще один ресурс, на который стоит обратить внимание: эта SIMD (SSE3 / SSE4) реализация 3D свертки Intel, который включает в себя как источник и презентация. Код для 16-битных целых чисел. Если вы не перейдете в GPU (например, CUFFT), это, вероятно, трудно получить быстрее, чем реализации Intel, которая также включает в себя Intel MKL. Ниже приведен пример трехмерной свертки (поплавок одинарной точности) эта страница документации MKL (ссылка исправлена, теперь отображается в https://stackoverflow.com/a/27074295/2778484).

15

Другие решения

Вы можете попробовать методы overlap-add и overlap-save. Они включают разбиение вашего входного сигнала на более мелкие куски, а затем с помощью любого из вышеуказанных методов.

FFT наиболее вероятен — и я могу ошибаться — самый быстрый метод, особенно если вы используете встроенные подпрограммы в MATLAB или библиотеку в C ++. Кроме того, неплохо было бы разбить входной сигнал на более мелкие куски.

2

у меня есть 2 способа рассчитать fastconv

и 2 лучше, чем 1

1- броненосец
Вы можете использовать библиотеку броненосца для вызова конвона с этим кодом

cx_vec signal(1024,fill::randn);
cx_vec code(300,fill::randn);
cx_vec ans = conv(signal,code);

2 — используйте fftw ans sigpack и библиотеку armadillo для вызова fast conv таким образом, вы должны инициализировать fft вашего кода в конструкторе

FastConvolution::FastConvolution(cx_vec inpCode)
{
filterCode = inpCode;
fft_w = NULL;
}cx_vec FastConvolution::filter(cx_vec inpData)
{
int length = inpData.size()+filterCode.size();
if((length & (length - 1)) == 0)
{

}
else
{
length = pow(2 , (int)log2(length) + 1);
}
if(length != fftCode.size())
initCode(length);

static cx_vec zeroPadedData;
if(length!= zeroPadedData.size())
{
zeroPadedData.resize(length);
}
zeroPadedData.fill(0);
zeroPadedData.subvec(0,inpData.size()-1) = inpData;cx_vec fftSignal = fft_w->fft_cx(zeroPadedData);
cx_vec mullAns = fftSignal % fftCode;
cx_vec ans = fft_w->ifft_cx(mullAns);
return ans.subvec(filterCode.size(),inpData.size()+filterCode.size()-1);
}

void FastConvolution::initCode(int length)
{
if(fft_w != NULL)
{
delete fft_w;
}
fft_w = new sp::FFTW(length,FFTW_ESTIMATE);
cx_vec conjCode(length,fill::zeros);
fftCode.resize(length);
for(int i = 0; i < filterCode.size();i++)
{
conjCode.at(i) = filterCode.at(filterCode.size() - i - 1);
}
conjCode = conj(conjCode);
fftCode = fft_w->fft_cx(conjCode);
}
0