python - 我应该如何将 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 函数，但我只需要对转换的一半进行 fft - 我的过滤器可以提前进行一次 fft，然后一遍又一遍地使用。

Answer 1

您不必返回np.float64 和hstack。您可以创建一个空的目标数组，形状与 sfft.rfft(Y) 和 sfft.rfft(X) 相同，然后创建一个 np.complex128 View 并用乘法结果填充此 View 。这将根据需要自动填充目标数组。
如果我重新学习你的例子:

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