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

user_1711271 0
Viewed 53

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.])
convolution fft numpy python
1 Answers

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.])
user_93321870
Related