Die Berechnung der Faltung zweier Funktionen mit FFT (FFTW)

Ich versuche zu beschleunigen, eine Berechnung für einen neuralen simulator mit Hilfe der FFT.

Ist die Gleichung:

(1) \sum(j=1 bis N) (w(i - j) * s_NMDA[j])

wo s_NMDA ist ein Vektor der Länge N und w ist definiert durch:

(2) w(j) = tanh[1/(2 * sigma * p)] * exp(-abs(j) /(sigma * p)]

wo sigma und p sind Konstanten.

(gibt es eine bessere Möglichkeit zum Rendern von Formeln auf stackoverflow?)

Die Berechnung muss durchgeführt werden, für die N Neuronen. Da (1) hängt nur von der absoluten Distanz abs(i - j), sollte es möglich sein, zu berechnen, diese mit Hilfe der FFT (convolution theorem).

Habe ich versucht, diese umsetzen mit FFTW aber die Ergebnisse passen nicht mit den erwarteten Ergebnissen. Ich habe noch nie verwendet FFTW vor und ich bin mir jetzt nicht sicher, ob meine Implementierung falsch ist, wenn meine Annahmen über das faltungs-theorem false sind.

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);

Hier ist ein Beispiel für die Ausgabe für 20 Neuronen:

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

Scheinen die Ergebnisse zu sein, fast richtig, aber der Fehler steigt mit der Anzahl der Neuronen. Auch scheinen die Ergebnisse um genauer zu sein für die Positionen (i), die sehr niedrig oder sehr hoch. Was ist denn hier Los?

Update: Wie vorgeschlagen von Oli Charlesworth, ich implementierte den Algorithmus in octave zu sehen, ob es eine Umsetzung oder Mathe-problem:

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)

Ergebnisse:

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

Sind die Ergebnisse sehr ähnlich. Daher muss es etwas falsch mit meinem Verständnis der convolution theorem /FFT.

StackOverflow hat keine direkte Unterstützung für Latex. Sie können jedoch verwenden Sie eine Website wie codecogs.com/latex/eqneditor.php, die dynamisch erzeugt Bilder.
Für den Kern des Problems, denke ich, müssen Sie zuerst, um festzustellen, ob Sie ein Mathe-problem oder ein problem FFTW. Ich würde vorschlagen, prototyping Ihren Algorithmus in so etwas wie Matlab/Octave. Alternativ könnten Sie versuchen, die Durchführung der Berechnung mit einem naiven DFT-Algorithmus (es ist etwa 7 Zeilen code). Wenn das funktioniert, dann wissen Sie, machen Sie etwas falsch mit FFTW!
Vielen Dank für Ihre Anregungen. Leider kann ich keine Bilder, weil ich nicht genug Ruf noch. Ich werde versuchen, das zu implementieren einen einfachen Prototyp in die Oktave und die Ergebnisse vergleichen.

InformationsquelleAutor nebw | 2012-06-24

biological-neural-network c++fft fftw math

8

To convolve 2 Signale mittels FFT im Allgemeinen werden Sie brauchen, um dies zu tun:
1. Fügen Sie so viele Nullen angehängt, um jedes signal als notwendig, damit seine Länge wird die kumulative Länge der original-Signale - 1 (das ist die Länge des Ergebnis der Faltung).
2. Wenn Ihr FFT-Bibliothek erfordert Eingabe Längen auf Vielfache von 2 hinzufügen, um jedes signal, wie viele Nullen wie nötig, um dieser Anforderung gerecht zu.
3. Berechnen Sie die DFT des Signals 1 (über FFT).
4. Berechnen Sie die DFT des Signals 2 (über FFT).
5. Multiplizieren Sie die beiden DFTs element-Weise. Es soll eine komplexe Multiplikation, btw.
6. Berechnen Sie die inverse DFT (über FFT) der multipliziert DFTs. Dass Sie Ihre convolution Ergebnis.
In deinem code sehe ich FFTW_FORWARD in allen 3 FFTs. Ich vermute, wenn das ist nicht das problem, es ist ein Teil davon. Die Letzte FFT sollte "rückwärts", nicht "vorwärts".

Außerdem denke ich, müssen Sie "+" in der 2. Ausdruck nicht "-":
```
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;
```
- Vielen Dank für Ihre Hilfe. Sie sind richtig über die zwei Fehler in meinem code. Sie waren nur ein Teil des Problems, aber. Es gab ein problem mit der Reihenfolge der Entfernungen array, aber ich kann keine Antwort auf meine eigene Frage für die nächsten fünf Stunden, da ich weniger als 10 Ruf, also ich denke, ich werde warten müssen, bevor ich es poste...
InformationsquelleAutor Alexey Frunze
1

Ich endlich das problem gelöst, vielen Dank an Alex und Oli Charlesworth für Eure Vorschläge!
```
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)
```
Ergebnisse:
```
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
```
Ich im Grunde vergessen, um die Entfernungen array in der richtigen Art und Weise. Wenn jemand interessiert ist, kann ich eine genauere Erklärung später.

Update: (Erklärung)

Hier ist, was mein Entfernungen vector (5 Neuronen) sah zunächst:
```
i =  1       2       3       4       5

| _w(0) | _w(1) | _w(2) | _w(3) | _w(4) |
```
Auf diesem Vektor, wandte ich ein kernel, z.B.:
```
|  0.1  |  0.1  |  0.0  |  0.2  |  0.3  |
```
Da war ich mit der zyklischen Faltung, das Ergebnis für das erste neuron _w(0) war:

0.0 * _w(2) + 0.1 * _w(1) + 0.1 * _w(0) + 0.1 * _w(1) + 0.0 * _w(2)

Aber das ist falsch, das Ergebnis sollte sein:

0.1 * _w(0) + 0.1 * _w(1) + 0.0 * _w(2) + 0.2 * _w(3) + 0.3 * _w(4)

Zu erreichen, musste ich "Spiegel" meine Strecken Vektor-und fügen Sie etwas Polsterung an den kernel:
```
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    |
```
Nun, habe ich, wenn ich die Faltung, die Ergebnisse für i = [5:9] sind genau die Ergebnisse, die ich suchte, also werde ich nur haben, um Sie zu verwerfen [1:4] und ich bin fertig 🙂
```
convolution = ifft(convolution_f)(n:end)
```
- So, Ihre Lösung, nicht um C (++) - code und FFTW nicht mehr? Ich denke, was ich versuche zu sagen, eine Antwort auf eine Frage auf StackOverflow sollte auch die Antwort auf die ursprüngliche Frage und nicht nur eine bloße "vielen Dank, es funktioniert jetzt", das ist das, was diese "Antwort" sieht. Sie können Ihre eigene Frage zu beantworten, das ist in Ordnung, aber es sollte trotzdem eine Antwort, nicht nur ein Gespräch. Für die, die wir haben Kommentare, in denen diese "Antwort" ist, glaube ich, gehört. Wenn Sie schon geholfen " Antworten, bitte vote und/oder akzeptieren, was auch immer angemessen ist. Das ist, wie SO funktioniert.
- Hey, wie bereits im update zu meiner Frage, das problem stellte sich heraus zu sein, ein Mathe-problem und nicht ein problem FFTW als ich ursprünglich dachte. Während Alex 's Antwort war wirklich hilfreich für mich, es didn' T lösen das eigentliche problem, das ich konfrontiert war, so wollte ich nicht akzeptieren, seine Antwort. Ich wollte upvote ihn, aber ich habe nicht genug Ruf zu tun. Der Grund, warum meine Antwort umfasst nur matlab-code und kein C++ - code ist, weil die Menschen können versuchen, den matlab-code für sich selbst, während der C++ - code in meine ursprüngliche Frage ist ein Teil eines größeren Programms und läuft nicht alleine.
- Bitte korrigieren Sie mich, wenn ich falsch Liege, ich bin neu in SO. Kann ich auch per post eine aktualisierte C++ version, aber ich weiß nicht wirklich den Punkt sehen, da, abgesehen von der zwei Probleme, die Alex erwähnt, es ist nichts falsch mit dem code, es funktioniert einfach nicht, was ich erwartet hatte, es zu tun zunächst.
- Der Punkt ist, wenn niemand die Antwort wirklich geholfen hat, oder Sie haben die beste Antwort, du bist frei, um Ihre eigene Frage zu beantworten (und auch akzeptieren). Aber ich glaube nicht, das ist hier der Fall, wie Sie wurden, am spitz 2 große Fehler in Ihrem C++ - Implementierung und der Bestätigung durch andere Mittel (was war das für eine Sprache?) dass der Algorithmus korrekt ist, und es ist nur die C++ - Implementierung, die ist falsch. In diesem Fall, IMHO, deine Antwort sollte gelöscht werden, und mir akzeptiert.
- Eigentlich ist das, was ich bestätigte mit meiner Matlab/Octave Umsetzung war, dass mein Algorithmus war nicht korrekt sind und der Fehler in meinem C++ - Implementierung waren nur ein Teil des Problems. Dehnte ich meine Antwort so, dass es erklärt, wie gründlich meine aktuelle Lösung. Ich denke, ich könnte noch markieren Sie Ihre Antwort als akzeptiert (da es noch half mir lösen meine C++ Implementierung), aber die eigentliche Lösung für mein problem ist in meiner Antwort.
- Ich sehe. Gute Erklärung. Gehen Sie wie Sie möchten. 🙂
- Bitte posten Sie Ihren C-code. Wenn du es in C, das wäre sehr hilfreich, da so weit ich erzählen kann, gibt es keine Beispiele, wie es zu tun in C. ich verstehe nicht, Ihren Widerstand. Sie sollte die post eine Lösung für Ihr problem, da Sie gestellt war, nicht von Ihr wahrgenommen-Lösung. Was Sie getan haben, ist yeild Menschen wie mich, die kommen später, die sind frustriert und haben zu bitten, eine separate Frage.
InformationsquelleAutor nebw

Schreibe einen Kommentar

Du musst angemeldet sein, um einen Kommentar abzugeben.