使用CUDA实现希尔伯特变换

Question

为了对一维数组执行Hilbert变换，需要执行以下步骤：

• 对数组进行FFT
• 将数组的一半加倍，另一半清零
• 对结果进行反FFT

我正在使用PyCuLib进行FFT。 我目前的代码：

def htransforms(data):
    N = data.shape[0]
    transforms       = nb.cuda.device_array_like(data)          # Allocates memory on GPU with size/dimensions of signal
    transforms.dtype = np.complex64                             # Change GPU array type to complex for FFT 

    pyculib.fft.fft(signal.astype(np.complex64), transforms)    # Do FFT on GPU

    transforms[1:N/2]      *= 2.0      # THIS STEP DOESN'T WORK
    transforms[N/2 + 1: N]  = 0+0j     # NEITHER DOES THIS ONE

    pyculib.fft.ifft_inplace(transforms)                        # Do IFFT on GPU: in place (same memory)    
    envelope_function      = transforms.copy_to_host()          # Copy results to host (computer) memory
    return abs(envelope_function)

我有一种感觉，这可能与Numba的CUDA接口本身有关...它是否允许像这样修改数组的单个元素（或数组切片）？我认为它可能会，因为变量transforms是一个numba.cuda.cudadrv.devicearray.DeviceNDArray，所以我想它可能具有一些与numpy的ndarray相同的操作。

简而言之，使用Numba的device_arrays，如何在切片上执行简单操作？我得到的错误是：

unsupported operand type(s) for *=: 'DeviceNDArray' and 'float'

- Yoshi

---

data 是什么类型？你对数组的 dtype 的操作看起来非常可疑。 - talonmies

更改数据类型不会改变数据本身。实际上，它什么也不做，也没有任何意义。双精度复数的数据类型为complex128，而单精度复数的数据类型为complex64。如果需要在不同类型之间进行转换，请使用astype。 - talonmies

很酷，我改了几行代码来考虑这个问题，但是我遇到的主要问题是“不支持的操作数类型 *=: 'DeviceNDArray' 和'float'”，所以我想知道是否有一种方法可以直接在GPU内存中修改此数组。 - Yoshi

由于它是一个“DeviceNDArray”，如果您想在GPU上完成所有操作，您可能需要编写一个CUDA内核来执行双倍和零位操作。在http://numba.pydata.org/numba-doc/0.20.0/cuda/kernels.html#absolute-positions上的示例似乎非常简单明了。 - AKX

或者，你可以在数组的切片上调用 cuBLAS.scal(2) 和 cuBLAS.scal(0)。http://pyculib.readthedocs.io/en/latest/cublas.html - AKX

Answer 1

我会使用pytorch来完成

您的函数使用pytorch（并且我删除了abs以返回复数）

def htransforms(data):
    N = data.shape[-1]
    # Allocates memory on GPU with size/dimensions of signal
    transforms       = torch.tensor(data).cuda()
    
    torch.fft.fft(transforms, axis=-1)
    transforms[:, 1:N//2]      *= 2.0      # THIS STEP DOESN'T WORK
    transforms[:, N//2 + 1: N]  = 0+0j     # NEITHER DOES THIS ONE

    # Do IFFT on GPU: in place (same memory)
    return torch.abs(torch.fft.ifft(transforms)).cpu() 

但是你的转换实际上与我在wikipedia上找到的不同

维基百科版本

def htransforms_wikipedia(data):
    N = data.shape[-1]
    # Allocates memory on GPU with size/dimensions of signal
    transforms       = torch.tensor(data).cuda()
    
    transforms = torch.fft.fft(transforms, axis=-1)
    transforms[:, 1:N//2]      *= -1j      # positive frequency
    transforms[:, (N+2)//2 + 1: N] *= +1j # negative frequency
    transforms[:,0] = 0; # DC signal
    if N % 2 == 0:
        transforms[:, N//2] = 0; # the (-1)**n term
    
    # Do IFFT on GPU: in place (same memory)
    return torch.fft.ifft(transforms).cpu() 

data = torch.zeros((1, 2**10))
data[:, 2**9] = 1;
tdata = htransforms(data).data;
plt.plot(tdata.real.T, '-')
plt.plot(tdata.imag.T, '-')
plt.xlim([500, 525])
plt.legend(['real', 'imaginary'])
plt.title('inpulse response of your version')

data = torch.zeros((1, 2**10))
data[:, 2**9] = 1;
tdata = htransforms_wikipedia(data).data;
plt.plot(tdata.real.T, '-');
plt.plot(tdata.imag.T, '-');
plt.xlim([500, 525])
plt.legend(['real', 'imaginary'])
plt.title('inpulse response of Wikipedia version')

你的版本的脉冲响应是1 + 1j * h[k]，其中h[k]是维基百科版本的脉冲响应。如果你正在处理真实数据，维基百科版本很好，因为你可以使用rfft和irfft生成一个线性版本。

def real_htransforms_wikipedia(data):
    N = data.shape[-1]
    # Allocates memory on GPU with size/dimensions of signal
    transforms       = torch.tensor(data).cuda()
    transforms = -1j * torch.fft.rfft(transforms, axis=-1)
    transforms[0] = 0;
    return torch.fft.irfft(transforms, axis=-1)

- Bob