Вычисление свертки двух функций с использованием FFT (FFTW)

Я пытаюсь ускорить вычисления для нейронного симулятора с использованием БПФ.

Уравнение:

(1) \ сумма (от j = 1 до N) (w (i-j) * s_NMDA [j])

где s_NMDA - вектор длины N, а w определяется как:

(2) w (j) = tanh [1 / (2 * сигма * p)] * exp (-abs (j) / (сигма * p)]

где сигма и р постоянные.

(Есть ли лучший способ визуализации уравнений на стеке потока?)

Расчет должен быть сделан для N нейронов. Поскольку (1) зависит только от абсолютного расстояния abs (i - j), его можно рассчитать с использованием БПФ (теорема свертки).

Я пытался реализовать это с использованием FFTW, но результаты не совпадают с ожидаемыми. Я никогда не использовал FFTW раньше, и теперь я не уверен, что моя реализация неверна, если мои предположения о теореме о свертке неверны.

void f_I_NMDA_FFT(
    const double     **states, // states[i][6] == s_NMDA[i]
    const unsigned int numNeurons)
{
    fftw_complex *distances, *sNMDAs, *convolution;
    fftw_complex *distances_f, *sNMDAs_f, *convolution_f;
    fftw_plan     p, pinv;
    const double scale = 1./numNeurons;

        distances = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * numNeurons);
        sNMDAs    = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * numNeurons);
        convolution = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * numNeurons);
        distances_f = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * numNeurons);
        sNMDAs_f    = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * numNeurons);
        convolution_f    = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * numNeurons);

        // fill input array for distances  
    for (unsigned int i = 0; i < numNeurons; ++i)
    {
        distances[i][0] = w(i);
        distances[i][1] = 0;
    }

        // fill input array for sNMDAs
    for (unsigned int i = 0; i < numNeurons; ++i)
    {
        sNMDAs[i][0] = states[i][6];
        sNMDAs[i][1] = 0;
    }

    p = fftw_plan_dft_1d(numNeurons,
                         distances,
                         distances_f,
                         FFTW_FORWARD,
                         FFTW_ESTIMATE);
    fftw_execute(p);

    p = fftw_plan_dft_1d(numNeurons,
                         sNMDAs,
                         sNMDAs_f,
                         FFTW_FORWARD,
                         FFTW_ESTIMATE);
    fftw_execute(p);

    // convolution in frequency domain
    for(unsigned int i = 0; i < numNeurons; ++i)
    {
        convolution_f[i][0] = (distances_f[i][0] * sNMDAs_f[i][0]
            - distances_f[i][1] * sNMDAs_f[i][1]) * scale;
        convolution_f[i][1] = (distances_f[i][0] * sNMDAs_f[i][1]
            - distances_f[i][1] * sNMDAs_f[i][0]) * scale;
    }

    pinv = fftw_plan_dft_1d(numNeurons,
                            convolution_f,
                            convolution,
                            FFTW_FORWARD,
                            FFTW_ESTIMATE);
    fftw_execute(pinv);

    // compute and compare with expected result
    for (unsigned int i = 0; i < numNeurons; ++i)
    {
            double expected = 0;

            for (int j = 0; j < numNeurons; ++j)
            {
                expected += w(i - j) * states[j][6];
            }
            printf("i=%d, FFT: r%f, i%f : Expected: %f\n", i, convolution[i][0], convolution[i][1], expected);
    }

    fftw_destroy_plan(p);
    fftw_destroy_plan(pinv);

    fftw_free(distances), fftw_free(sNMDAs), fftw_free(convolution);
    fftw_free(distances_f), fftw_free(sNMDAs_f), fftw_free(convolution_f);

Вот пример вывода для 20 нейронов:

i=0, FFT: r0.042309, i0.000000 : Expected: 0.041504
i=1, FFT: r0.042389, i0.000000 : Expected: 0.042639
i=2, FFT: r0.042466, i0.000000 : Expected: 0.043633
i=3, FFT: r0.042543, i0.000000 : Expected: 0.044487
i=4, FFT: r0.041940, i0.000000 : Expected: 0.045203
i=5, FFT: r0.041334, i0.000000 : Expected: 0.045963
i=6, FFT: r0.041405, i0.000000 : Expected: 0.046585
i=7, FFT: r0.041472, i0.000000 : Expected: 0.047070
i=8, FFT: r0.041537, i0.000000 : Expected: 0.047419
i=9, FFT: r0.041600, i0.000000 : Expected: 0.047631
i=10, FFT: r0.041660, i0.000000 : Expected: 0.047708
i=11, FFT: r0.041717, i0.000000 : Expected: 0.047649
i=12, FFT: r0.041773, i0.000000 : Expected: 0.047454
i=13, FFT: r0.041826, i0.000000 : Expected: 0.047123
i=14, FFT: r0.041877, i0.000000 : Expected: 0.046656
i=15, FFT: r0.041926, i0.000000 : Expected: 0.046052
i=16, FFT: r0.041294, i0.000000 : Expected: 0.045310
i=17, FFT: r0.042059, i0.000000 : Expected: 0.044430
i=18, FFT: r0.042144, i0.000000 : Expected: 0.043412
i=19, FFT: r0.042228, i0.000000 : Expected: 0.042253

Результаты кажутся почти правильными, но ошибка увеличивается с числом нейронов. Кроме того, результаты кажутся более точными для позиций (i), которые являются очень низкими или очень высокими. Что здесь происходит?

UpdateКак предложил Оли Чарльзуорт, я реализовал алгоритм в октаве, чтобы увидеть, является ли это реализацией или математической проблемой:

input = [0.186775; 0.186775; 0.186775; 0.186775; 0.186775; 0; 0.186775; 0.186775; 0.186775; 0.186775];

function ret = _w(i)
  ret = tanh(1 / (2* 1 * 32)) * exp(-abs(i) / (1 * 32));
end

for i = linspace(1, 10, 10)
  expected = 0;
  for j = linspace(1, 10, 10)
    expected += _w(i-j) * input(j);
  end
  expected
end

distances = _w(transpose(linspace(0, 9, 10)));

input_f = fft(input);
distances_f = fft(distances);

convolution_f = input_f .* distances_f;

convolution = ifft(convolution_f)

Результаты:

expected =  0.022959
expected =  0.023506
expected =  0.023893
expected =  0.024121
expected =  0.024190
expected =  0.024100
expected =  0.024034
expected =  0.023808
expected =  0.023424
expected =  0.022880
convolution =

   0.022959
   0.023036
   0.023111
   0.023183
   0.023253
   0.022537
   0.022627
   0.022714
   0.022798
   0.022880

Результаты очень похожи. Следовательно, должно быть что-то не так с моим пониманием теоремы свертки / БПФ.

 nebw24 июн. 2012 г., 16:48
Спасибо за ваши предложения. К сожалению, я не могу включить изображения, потому что у меня еще недостаточно репутации. Я попытаюсь реализовать простой прототип в октаве и сравнить результаты.
 Oliver Charlesworth24 июн. 2012 г., 16:32
Что касается основной проблемы, я думаю, что вам прежде всего необходимо установить, рассматриваете ли вы математическую проблему или проблему FFTW. Я хотел бы предложить прототип вашего алгоритма в нечто вроде Matlab / Octave. В качестве альтернативы вы могли бы попытаться выполнить вычисление с помощью простого алгоритма DFT (он содержит около 7 строк кода). Если это работает, то вы знаете, что вы делаете что-то не так с FFTW!
 Oliver Charlesworth24 июн. 2012 г., 16:29
StackOverflow не имеет прямой поддержки латекса. Тем не менее, вы можете использовать сайт, какcodecogs.com/latex/eqneditor.php, который динамически генерирует изображения.

Ответы на вопрос(2)

Решение Вопроса

большое спасибо Алексу и Оли Чарльзуорту за ваши предложения!

function ret = _w(i)
  ret = tanh(1 / (2* 1 * 32)) * exp(-abs(i) / (1 * 32));
end

_input = [0.186775; 0.186775; 0.186775; 0.186775; 0.186775; 0; 0.186775; 0.186775; 0.186775; 0.186775];
n = size(_input)(1);

input = _input;

for i = linspace(1, n, n)
  expected = 0;
  for j = linspace(1, n, n)
    expected += _w(i-j) * input(j);
  end
  expected
end

input = vertcat(_input, zeros((2*n)-n-1,1));

distances = _w(transpose(linspace(0, (2*n)-n-1, n)));
distances = vertcat(flipud(distances), distances(2:end));

input_f = fft(input);
distances_f = fft(distances);

convolution_f = input_f .* distances_f;

convolution = ifft(convolution_f)(n:end)

Результаты:

expected =  0.022959
expected =  0.023506
expected =  0.023893
expected =  0.024121
expected =  0.024190
expected =  0.024100
expected =  0.024034
expected =  0.023808
expected =  0.023424
expected =  0.022880
convolution =

   0.022959
   0.023506
   0.023893
   0.024121
   0.024190
   0.024100
   0.024034
   0.023808
   0.023424
   0.022880

Я в основном забыл правильно упорядочить массив расстояний. Если кому-то интересно, я могу предоставить более подробное объяснение позже.

Update: (Explanation)

Вот то, что мой вектор расстояний (для 5 нейронов) выглядел изначально так:

i =  1       2       3       4       5

| _w(0) | _w(1) | _w(2) | _w(3) | _w(4) |

На этот вектор я применил ядро, например:

|  0.1  |  0.1  |  0.0  |  0.2  |  0.3  |

Поскольку я использовал циклическую свертку, результат для первого нейрона _w (0) был:

0,0 * _w (2) + 0,1 * _w (1) + 0,1 * _w (0) + 0,1 * _w (1) + 0,0 * _w (2)

Но это неверно, результат должен быть:

0,1 * _w (0) + 0,1 * _w (1) + 0,0 * _w (2) + 0,2 * _w (3) + 0,3 * _w (4)

Чтобы достичь этого, мне нужно было «отразить» мой вектор расстояний и добавьте немного дополнения к ядру:

input = vertcat(_input, zeros((2*n)-n-1,1));
distances = _w(transpose(linspace(0, (2*n)-n-1, n)));
distances = _w(transpose(linspace(0, (2*n)-n-1, n)));

i =  1       2        3       4       5       6       7       8       9

| _w(4) | _w(3)  | _w(2) | _w(1) | _w(0) | _w(1) | _w(2) | _w(3) | _w(4) |

|  0.1  |  0.1   |  0.0  |  0.2  |  0.3  |  0    |  0    |  0    |  0    |

Теперь, если я применю свертку, то результаты для i = [5: 9] будут именно теми результатами, которые я искал, поэтому мне нужно только отказаться от [1: 4], и я сделал :)

convolution = ifft(convolution_f)(n:end)
 nebw25 июн. 2012 г., 11:45
На самом деле, то, что я подтвердил с моей реализацией Matlab / Octave, это то, что мой алгоритмnot Правильно, и что ошибки в моей реализации C ++ были только частью проблемы. Я расширил свой ответ, чтобы он более подробно объяснил мое реальное решение. Я думаю, что я все еще мог бы пометить ваш ответ как принятый (так как он все еще помог мне исправить мою реализацию C ++), но фактическое решение моей проблемы находится в моем ответе.
 nebw25 июн. 2012 г., 11:09
Пожалуйста, исправьте меня, если я ошибаюсь, я новичок в SO. Я также могу опубликовать обновленную версию C ++, но на самом деле я не вижу в этом смысла, поскольку, помимо двух упомянутых Алексом проблем, в коде нет ничего плохого, он просто не выполняет то, что ожидал. первоначально.
 25 июн. 2012 г., 00:47
Итак, ваше решение больше не включает в себя код C (++) и FFTW? Я предполагаю, что то, что я пытаюсь сказать, ответ на вопрос в StackOverflow должен отвечать на исходный вопрос, а не просто быть «спасибо, он работает сейчас», на что это «ответ». выглядит. Вы можете ответить на свой собственный вопрос, это нормально, но, тем не менее, это должен быть ответ, а не просто разговор. Для этого у нас есть комментарии, к которым этот «ответ», я полагаю, принадлежит. Если вам помогли другие? Ответы, пожалуйста, оцените их и / или примите, что угодно. Вот как это работает.
 25 июн. 2012 г., 11:14
Дело в том, что если никто не помог вам, или у вас есть лучший ответ, вы можете ответить на свой вопрос (и даже принять его). Но я не думаю, что это имеет место в данном случае, поскольку вы указали на две основные ошибки в вашей реализации C ++ и другими средствами (что это был за язык?) Вы подтвердили, что алгоритм правильный и это только реализация C ++. это неправильно. В этом случае, ИМХО, ваш ответ должен быть удален, а мой принят.
 nebw25 июн. 2012 г., 11:00
Эй, как указано в обновлении моего вопроса, проблема оказалась математической, а не FFTW-проблемой, как я изначально думал. Хотя ответ Алекса был действительно полезным для меня, он не решил реальную проблему, с которой я столкнулся, поэтому я не принял его ответ. Я хотел бы поддержать его, но у меня недостаточно репутации, чтобы сделать это. Причина, по которой мой ответ включает только код matlab и никакой код C ++, состоит в том, что люди могут сами попробовать код matlab, в то время как код C ++ в моем исходном вопросе является частью более крупной программы и не запускается на ее основе.

как правило, нужно сделать следующее:

Add as many zeroes to every signal as necessary so its length becomes the cumulative length of the original signals - 1 (that's the length of the result of the convolution). If your FFT library requires input lengths to be powers of 2, add to every signal as many zeroes as necessary to meet that requirement too. Calculate the DFT of signal 1 (via FFT). Calculate the DFT of signal 2 (via FFT). Multiply the two DFTs element-wise. It should be a complex multiplication, btw. Calculate the inverse DFT (via FFT) of the multiplied DFTs. That'll be your convolution result.

В вашем коде я вижуFFTW_FORWARD во всех 3 БПФ. Я догадываюсь, если это не проблема, это ее часть. Последнее БПФ должно быть «назад», а не «вперед».

Кроме того, я думаю, что вам нужно & quot; + & quot; во втором выражении, а не "-":

convolution_f[i][0] = (distances_f[i][0] * sNMDAs_f[i][0]
            - distances_f[i][1] * sNMDAs_f[i][1]) * scale;

convolution_f[i][1] = (distances_f[i][0] * sNMDAs_f[i][1]
            - distances_f[i][1] * sNMDAs_f[i][0]) * scale;
 nebw24 июн. 2012 г., 18:58
Спасибо за вашу помощь. Вы правы в отношении двух ошибок в моем коде. Хотя они были только частью проблемы. Была другая проблема с порядком массива расстояний, но я не могу ответить на свой вопрос в течение следующих пяти часов, потому что у меня менее 10 репутаций, поэтому я думаю, что мне придется подождать, прежде чем я смогу опубликовать его. ,

Ваш ответ на вопрос