Databases
StackOverflow2013
Postsdbo
POConvolution computations in Numpy/Scipy
 

Note that there are some explanatory texts on larger screens.

plurals
  1. POConvolution computations in Numpy/Scipy
    primarykey
    data
    text
    <p>Profiling some computational work I'm doing showed me that one bottleneck in my program was a function that basically did this (<code>np</code> is <code>numpy</code>, <code>sp</code> is <code>scipy</code>):</p> <pre><code>def mix1(signal1, signal2): spec1 = np.fft.fft(signal1, axis=1) spec2 = np.fft.fft(signal2, axis=1) return np.fft.ifft(spec1*spec2, axis=1) </code></pre> <p>Both signals have shape <code>(C, N)</code> where <code>C</code> is the number of sets of data (usually less than 20) and <code>N</code> is the number of samples in each set (around 5000). The computation for each set (row) is completely independent of any other set.</p> <p>I figured that this was just a simple convolution, so I tried to replace it with:</p> <pre><code>def mix2(signal1, signal2): outputs = np.empty_like(signal1) for idx, row in enumerate(outputs): outputs[idx] = sp.signal.convolve(signal1[idx], signal2[idx], mode='same') return outputs </code></pre> <p>...just to see if I got the same results. But I didn't, and my questions are:</p> <ol> <li>Why not?</li> <li>Is there a better way to compute the equivalent of <code>mix1()</code>?</li> </ol> <p>(I realise that <code>mix2</code> probably wouldn't have been faster as-is, but it might have been a good starting point for parallelisation.)</p> <p>Here's the full script I used to quickly check this:</p> <pre><code>import numpy as np import scipy as sp import scipy.signal N = 4680 C = 6 def mix1(signal1, signal2): spec1 = np.fft.fft(signal1, axis=1) spec2 = np.fft.fft(signal2, axis=1) return np.fft.ifft(spec1*spec2, axis=1) def mix2(signal1, signal2): outputs = np.empty_like(signal1) for idx, row in enumerate(outputs): outputs[idx] = sp.signal.convolve(signal1[idx], signal2[idx], mode='same') return outputs def test(num, chans): sig1 = np.random.randn(chans, num) sig2 = np.random.randn(chans, num) res1 = mix1(sig1, sig2) res2 = mix2(sig1, sig2) np.testing.assert_almost_equal(res1, res2) if __name__ == "__main__": np.random.seed(0x1234ABCD) test(N, C) </code></pre>
    singulars
    1. PO
      singulars
      1. PTAnswer
    1. This table or related slice is empty.
    1. PTQuestion
    1. This table or related slice is empty.
    1. USdetly
    plurals
    1. This table or related slice is empty.
    1. PL
      singulars
      1. LTLinked
    2. PL
      singulars
      1. LTLinked
    1. This table or related slice is empty.
    1. PO
      singulars
      1. PTAnswer
    2. PO
      singulars
      1. PTAnswer
    3. PO
      singulars
      1. PTAnswer
    1. VO
      singulars
      1. POConvolution computations in Numpy/Scipy
      1. This table or related slice is empty.
      1. VTUpMod
    2. VO
      singulars
      1. POConvolution computations in Numpy/Scipy
      1. This table or related slice is empty.
      1. VTUpMod
    3. VO
      singulars
      1. POConvolution computations in Numpy/Scipy
      1. This table or related slice is empty.
      1. VTUpMod
    1. COI don't know for sure but I suspect the difference may be due to the fact that FFT/IFFT are "circular", ie they make up for the fact that the data doesn't go on forever by pretending the signal is periodic, while convolve() fills in with zeros. There should be a way to due a circular convolution -- not sure how, though.
      singulars
      1. POConvolution computations in Numpy/Scipy
      1. USOwen
    2. CO@Owen - That must be it — `mix1` is doing the equivalent of a [circular convolution](http://en.wikipedia.org/wiki/Circular_convolution). I completely forgot about that distinction.
      singulars
      1. POConvolution computations in Numpy/Scipy
      1. USdetly
    3. COThere is also another difference because of normalization -- FFT normalizes its output by dividing by I think the length of the array (not all FFT implementations due this, but Numpy's does), whereas convolve doesn't due that because it thinks you actually want to add stuff up, not average it. Edit -- not sure about this.
      singulars
      1. POConvolution computations in Numpy/Scipy
      1. USOwen
 

Querying!

 
Guidance
A row detail

Detail views are divided into sections. All the information in the data section comes from columns in the selected row. The other sections display data from other, related rows.

Related data can be related in a to-one or a to-many fashion. Captions of data related in a to-many fashion link to a list view showing a filtered view of the table.

Try moving around until you find a non-empty to-many entry and click on the label to get to one. You can move back to the root by clicking on the database name in the header.

SQuiL has stopped working due to an internal error.

If you are curious you may find further information in the browser console, which is accessible through the devtools (F12).

Reload