我應該如何將scipy.fftpack輸出向量放在一起？

Question

scipy.fftpack.rfft函數返回DFT作爲浮點數的矢量，在實數和複數部分之間交替。這意味着一起乘以DFT（卷積），我將不得不做複雜的乘法「手動」，這似乎很棘手。這一定是人們經常做的事情 - 我認爲/希望有一個簡單的技巧來有效地做到這一點，我沒有發現？

基本上我想要使這兩種方法給出了相同的答案來解決這個代碼：

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

注：我知道fftpack包含一個convolve函數，但我只需要轉換一半的變換 - 我的過濾器可以提前一次，然後一遍又一遍地使用。

Answer 1

你可以把你的回報陣列的一個切片的視圖，例如：

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

您將需要以不同方式處理偶數或奇數大小的FFT：

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

Answer 2

你不必須翻轉回np.float64和hstack。您可以創建一個空目標數組，其形狀與sfft.rfft(Y)和sfft.rfft(X)相同，然後創建一個np.complex128視圖並填充此視圖以及乘法結果。這會自動填充目標數組。
如果我再拿你的例子：

import numpy as np 
import scipy.fftpack as sfft 

X = np.random.normal(size = 2000) 
Y = np.random.normal(size = 2000) 
Xf = np.fft.rfft(X) 
Xf_cpx = Xf[1:-1].view(np.complex128) 
Yf = np.fft.rfft(Y) 
Yf_cpx = Yf[1:-1].view(np.complex128) 

Zf = np.empty(X.shape) 
Zf_cpx = Zf[1:-1].view(np.complex128) 

Zf[0] = Xf[0]*Yf[0] 

# the [...] is important to use the view as a reference to Zf and not overwrite it 
Zf_cpx[...] = Xf_cpx * Yf_cpx 

Zf[-1] = Xf[-1]*Yf[-1] 

Z = sfft.irfft.irfft(Zf) 

，這就是它！ 你可以使用簡單的if語句，如果你希望你的代碼更一般，並處理奇數的長度，如Jaime的答案中所解釋的。 這是一個功能，你想要什麼：

def rfft_mult(a,b): 
    """Multiplies two outputs of scipy.fftpack.rfft""" 
    assert a.shape == b.shape 
    c = np.empty(a.shape) 
    c[...,0] = a[...,0]*b[...,0] 
    # To comply with the rfft support of multi dimensional arrays 
    ar = a.reshape(-1,a.shape[-1]) 
    br = b.reshape(-1,b.shape[-1]) 
    cr = c.reshape(-1,c.shape[-1]) 
    # Note that we cannot use ellipses to achieve that because of 
    # the way `view` work. If there are many dimensions, one should 
    # consider to manually perform the complex multiplication with slices. 
    if c.shape[-1] & 0x1: # if odd 
     for i in range(len(ar)): 
      ac = ar[i,1:].view(np.complex128) 
      bc = br[i,1:].view(np.complex128) 
      cc = cr[i,1:].view(np.complex128) 
      cc[...] = ac*bc 
    else: 
     for i in range(len(ar)): 
      ac = ar[i,1:-1].view(np.complex128) 
      bc = br[i,1:-1].view(np.complex128) 
      cc = cr[i,1:-1].view(np.complex128) 
      cc[...] = ac*bc 
     c[...,-1] = a[...,-1]*b[...,-1] 
    return c