Corone - 1 month ago 10x
Python Question

# How should I multiply scipy.fftpack output vectors together?

The

`scipy.fftpack.rfft`
function returns the DFT as a vector of floats, alternating between the real and complex part. This means to multiply to DFTs together (for convolution) I will have to do the complex multiplication "manually" which seems quite tricky. This must be something people do often - I presume/hope there is a simple trick to do this efficiently that I haven't spotted?

Basically I want to fix this code so that both methods give the same answer:

``````import numpy as np
import scipy.fftpack as sfft

X = np.random.normal(size = 2000)
Y = np.random.normal(size = 2000)
NZ = np.fft.irfft(np.fft.rfft(Y) * np.fft.rfft(X))
SZ = sfft.irfft(sfft.rfft(Y) * sfft.rfft(X))    # This multiplication is wrong

NZ
array([-43.23961083,  53.62608086,  17.92013729, ..., -16.57605207,
8.19605764,   5.23929023])
SZ
array([-19.90115323,  16.98680347,  -8.16608202, ..., -47.01643274,
-3.50572376,  58.1961597 ])
``````

N.B. I am aware that fftpack contains a
`convolve`
function, but I only need to fft one half of the transform - my filter can be fft'd once in advance and then used over and over again.

You can take a view of a slice of your return array, e.g.:

``````>>> scipy.fftpack.fft(np.arange(8))
array([ 28.+0.j        ,  -4.+9.65685425j,  -4.+4.j        ,
-4.+1.65685425j,  -4.+0.j        ,  -4.-1.65685425j,
-4.-4.j        ,  -4.-9.65685425j])
>>> a = scipy.fftpack.rfft(np.arange(8))
>>> a
array([ 28.        ,  -4.        ,   9.65685425,  -4.        ,
4.        ,  -4.        ,   1.65685425,  -4.        ])
>>> a.dtype
dtype('float64')
>>> a[1:-1].view(np.complex128) # First and last entries are real
array([-4.+9.65685425j, -4.+4.j        , -4.+1.65685425j])
``````

You will need to handle even or odd sized FFTs differently:

``````>>> scipy.fftpack.fft(np.arange(7))
array([ 21.0+0.j        ,  -3.5+7.26782489j,  -3.5+2.79115686j,
-3.5+0.79885216j,  -3.5-0.79885216j,  -3.5-2.79115686j,
-3.5-7.26782489j])
>>> a = scipy.fftpack.rfft(np.arange(7))
>>> a
array([ 21.        ,  -3.5       ,   7.26782489,  -3.5       ,
2.79115686,  -3.5       ,   0.79885216])
>>> a.dtype
dtype('float64')
>>> a[1:].view(np.complex128)
array([-3.5+7.26782489j, -3.5+2.79115686j, -3.5+0.79885216j])
``````