python - 使用 CUDA 进行希尔伯特变换
问题描述
为了对一维数组进行希尔伯特变换,必须:
- FFT 数组
- 将数组的一半加倍,另一半为零
- 逆FFT结果
我正在使用 PyCuLib 进行 FFTing。到目前为止我的代码
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
是 a numba.cuda.cudadrv.devicearray.DeviceNDArray
,所以我想它可能有一些与 numpy's 相同的操作ndarray
。
简而言之,使用 Numbadevice_arrays
对切片进行简单操作的最简单方法是什么?我得到的错误是
*= 不支持的操作数类型:“DeviceNDArray”和“float”
解决方案
我会使用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()
但是您的转换实际上与我在维基百科上发现的不同
维基百科版本
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)
推荐阅读
- c# - 在 C# 中强制计算四舍五入
- php - 如何从 php 数组中获取数据并使用 jquery 将其添加到页面标题属性中?
- c++ - C++ 中的 Json:将数字解析为字符串以避免浮点不准确
- java - 如何查看是否使用 Java 中的 PDFBox 检查了单选按钮或复选框?
- java - 在 Spring 中使用配置创建组件
- c# - 正则表达式向前看
- r - 模拟活跃权重在 [-0.02,0.02] 之间且绝对权重总和为 1 的投资组合
- javascript - Div 显示多次且按钮仍处于禁用状态
- javascript - 清单无法识别 jQuery
- c++ - 将指向成员函数的指针传递给模板