如何将scipy.fftpack输出向量相乘？ [英] How should I multiply scipy.fftpack output vectors together?

问题描述

scipy.fftpack.rfft 函数将DFT作为浮点向量返回，在实部和复杂部分之间交替。这意味着要一起乘以DFT（用于卷积），我将不得不手动执行复杂的乘法，这似乎很棘手。这一定是人们经常做的事情-我想/希望有一个简单的技巧可以有效地做到这一点，而我没有发现？

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

NB我知道fftpack包含一个 convolve 函数，但是我只需要fft转换的一半-我的过滤器可以提前fft'd一次，然后用于

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.

推荐答案

您不需要，您不必翻转回 np.float64 和 hstack 。您可以创建一个空的目标数组，其形状与 sfft.rfft（Y）和 sfft.rfft（X），然后为其创建一个 np.complex128 视图，并用乘法结果填充该视图。

如果我重新举一个例子：

You don't have to flip back to np.float64 and hstack. You can create an empty destination array, the same shape as sfft.rfft(Y) and sfft.rfft(X), then create a np.complex128 view of it and fill this view with the result of the multiplication. This will automatically fill the destination array as wanted.
If I retake your example :

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)

就是这样！
如果您希望代码更通用，并且可以按照Jaime的答案中所述处理奇数长度，则可以使用简单的if语句。
这是一个执行您想要的功能的函数：

and that's it! You can use a simple if statement if you want your code to be more general and handle odd lengths as explained in Jaime's answer. Here is a function that does what you want:

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

这篇关于如何将scipy.fftpack输出向量相乘？的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！