我使用 DHT 计算卷积的代码有什么问题？

Question

import numpy as np

def dht(x:np.array):    
    X = np.fft.fft(x)
    X = np.real(X) - np.imag(X)
    return X

def idht(X:np.array):
    n = len(X)
    X = dht(X)
    x = 1/n * X
    return x

根据 here, here(eq.4.7.22-4.7.25) 和许多其他地方说明的定理，可以使用 DHT 计算循环卷积。我实现了相应的算法并将其与使用 FFT 的循环卷积计算进行了比较，但结果并不接近。为什么？

def conv(x:np.array, y:np.array):
    X = dht(x)
    Y = dht(y)
    Xflip = np.flip(X)
    Yflip = np.flip(Y)
    Yplus = Y + Yflip
    Yminus = Y - Yflip
    Z = 0.5 * (X * Yplus + Xflip * Yminus)
    z = idht(Z)
    return z


def test_conv():
    x = np.ones((5, ))
    y = np.copy(x)
    z = conv(x, y)
    z1 = np.real(np.fft.ifft(np.fft.fft(x)*np.fft.fft(y)))
    np.testing.assert_allclose(z, z1, err_msg="test_convolution() failed")


if (__name__=='__main__'):
    test_conv()
    print("test_conv passed")

输出：

Traceback (most recent call last):
  File "ronf.py", line 35, in <module>
    test_conv()
  File "ronf.py", line 31, in test_conv
    np.testing.assert_allclose(z, z1, err_msg="test_convolution() failed")
  File "/home/andrea/vscode_venv/lib/python3.8/site-packages/numpy/testing/_private/utils.py", line 1528, in assert_allclose
    assert_array_compare(compare, actual, desired, err_msg=str(err_msg),
  File "/home/andrea/vscode_venv/lib/python3.8/site-packages/numpy/testing/_private/utils.py", line 842, in assert_array_compare
    raise AssertionError(msg)
AssertionError: 
Not equal to tolerance rtol=1e-07, atol=0
test_convolution() failed
Mismatched elements: 5 / 5 (100%)
Max absolute difference: 5.65018378
Max relative difference: 1.13003676
 x: array([ 0.      ,  4.105099,  5.992006,  3.053079, -0.650184])
 y: array([5., 5., 5., 5., 5.])

Answer 1

问题出在翻转部分。举例说明：

For a given array X_k, we need X_{-k} i.e., X_{N-k}.

N := 5
X_k := [0, 1, 2, 3, 4]

np.flip does this:
[4, 3, 2, 1, 0]

However, this is not X_{N-k}! We need:
X_{N-k} = [0, 4, 3, 2, 1]

since:
for k = 0 => X[5-0] = X[0] = 0
for k = 1 => X[5-1] = X[4] = 4
...                         ...
for k = 4 => X[5-4] = X[1] = 1

也就是说，我们假设周期为N s.t。 X[0] = X[N] 等等，但是 np.flip 的结果不符合那个。简而言之，第一个元素应该保持原样，其他元素应该翻转。

一种修复方法是 np.roll 在翻转后将其旋转 1 个位置：

def conv(x:np.array, y:np.array):
    X = dht(x)
    Y = dht(y)
    Xflip = np.roll(np.flip(X), shift=1)    # change is here
    Yflip = np.roll(np.flip(Y), shift=1)    # and here only
    Yplus = Y + Yflip
    Yminus = Y - Yflip
    Z = 0.5 * (X * Yplus + Xflip * Yminus)
    z = idht(Z)
    return z

然后我得到：

>>> a = np.ones((5,))
>>> conv(a, a)
array([5., 5., 5., 5., 5.])

我使用 DHT 计算卷积的代码有什么问题？

What is wrong in my code using DHT to compute convolution?

python

numpy

fft

convolution