15

Profiling some computational work I'm doing showed me that one bottleneck in my program was a function that basically did this (np is numpy, sp is scipy):

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)

Both signals have shape (C, N) where C is the number of sets of data (usually less than 20) and N is the number of samples in each set (around 5000). The computation for each set (row) is completely independent of any other set.

I figured that this was just a simple convolution, so I tried to replace it with:

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

...just to see if I got the same results. But I didn't, and my questions are:

  1. Why not?
  2. Is there a better way to compute the equivalent of mix1()?

(I realise that mix2 probably wouldn't have been faster as-is, but it might have been a good starting point for parallelisation.)

Here's the full script I used to quickly check this:

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)
Improve this question
5
  • I 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. Commented Jul 28, 2011 at 7:09
  • @Owen - That must be it — mix1 is doing the equivalent of a circular convolution. I completely forgot about that distinction. Commented Jul 28, 2011 at 7:11
  • There 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. Commented Jul 28, 2011 at 7:12
  • Also if you want to fix your bottleneck it'll help a lot to force your time axis to be a power of 2. Commented Jul 28, 2011 at 7:15
  • @Owen - huh, I thought that internally the FFT algorithm padded it up to the next power of 2. I'll check that. I'm not sure that Numpy's FFT does any normalisation — I remember checking something using Parseval's theorem, and I needed a factor of 1/N² in there. Not sure either. Commented Jul 28, 2011 at 7:19

3 Answers 3

Reset to default
10

So I tested this out and can now confirm a few things:

1) numpy.convolve is not circular, which is what the fft code is giving you:

2) FFT does not internally pad to a power of 2. Compare the vastly different speeds of the following operations:

x1 = np.random.uniform(size=2**17-1)
x2 = np.random.uniform(size=2**17)

np.fft.fft(x1)
np.fft.fft(x2)

3) Normalization is not a difference -- if you do a naive circular convolution by adding up a(k)*b(i-k), you will get the result of the FFT code.

The thing is padding to a power of 2 is going to change the answer. I've heard tales that there are ways to deal with this by cleverly using prime factors of the length (mentioned but not coded in Numerical Recipes) but I've never seen people actually do that.

Improve this answer
Sign up to request clarification or add additional context in comments.

8 Comments

Do you know, off the top of your head, if there's a "next power of 2" function in Numpy (like MATLAB's nextpow2)?
Haha no, I don't know. I've wanted one often as well :).
@detly If you expand both the orignal signals to 2^n and then convolve, you should get a real signal out. That also means you don't do any extra padding in the frequency domain, because that will break the symmetry.
@eryksun Well, that might be better, but it's still changing it. Say I have ABCDEF, and I want to extend it to 8 so it becomes ABCDEFAB. But that signal now represents ABCDEFABABCDEFAB... whereas the original was ABCDEFABCFED (the repeated AB is different).
Next power of two is very easy to compute: 2**ceil(log2(x))
|
2

scipy.signal.fftconvolve does convolve by FFT, it's python code. You can study the source code, and correct you mix1 function.

Improve this answer

1 Comment

The mix1 function is the one that works for me. I needed a mix2 function that replicated its behaviour.
1

As mentioned before, the scipy.signal.convolve function does not perform a circular convolution. If you want a circular convolution performed in realspace (in contrast to using fft's) I suggest using the scipy.ndimage.convolve function. It has a mode parameter which can be set to 'wrap' making it a circular convolution. 

for idx, row in enumerate(outputs):
    outputs[idx] = sp.ndimage.convolve(signal1[idx], signal2[idx], mode='wrap')
Improve this answer

1 Comment

@_ABDreverhaven I have tried what you said, but I am still getting different outputs

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.