Corone - 10 months ago 49

Python Question

The

`scipy.fftpack.rfft`

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`

Answer Source

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])
```