Самый быстрый метод расчета свертки
Кто-нибудь знает о самом быстром способе вычисления свертки? К сожалению, матрица, с которой я имею дело, очень большая (500x500x200), и если я использую convn в MATLAB это занимает много времени (мне приходится повторять этот расчет во вложенном цикле). Итак, я использовал свертку с FFT и теперь она быстрее. Но я все еще ищу более быстрый метод. Любая идея?
Решение
Если ваше ядро отделимо, наибольшее увеличение скорости будет реализовано при выполнении нескольких последовательных одномерных сверток.
Стив Эддинс из 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).
Другие решения
Вы можете попробовать методы overlap-add и overlap-save. Они включают разбиение вашего входного сигнала на более мелкие куски, а затем с помощью любого из вышеуказанных методов.
FFT наиболее вероятен — и я могу ошибаться — самый быстрый метод, особенно если вы используете встроенные подпрограммы в MATLAB или библиотеку в C ++. Кроме того, неплохо было бы разбить входной сигнал на более мелкие куски.
у меня есть 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);
}