The Discrete Hartley Transform can be computed as
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
Circular convolution can be computed using the DHT according to a theorem illustrated here, here(eq.4.7.22-4.7.25) and in many other places. I implemented the corresponding algorithm and compared it with circular convolution compute using the FFT, but the results are not even close. Why?
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")
Output:
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.])
The issue is about the flipping part. To show with an example:
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
That is, we assume periodicity of N s.t. X[0] = X[N] and so on, but np.flip's result doesn't obey that. In short, first element should stay as is, others should flip.
One way to fix is to np.roll after your flips to rotate it 1 position:
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
Then I get:
>>> a = np.ones((5,))
>>> conv(a, a)
array([5., 5., 5., 5., 5.])