# How to compute convolution using the Discrete Hartley Transform

It's easy to compute the Discrete Hartley Transform of a 1D signal:

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


The code for dht is based on this answer, while idht uses the involutory property of the DHT (see here) .

However, when I try to compute the circular convolution of two sequences using the well-known property of the DHT shown here, I get absurd results. What am I doing wrong? Here is the code...

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")


The code output is below: x is the circular convolution as computed by my code, and y is the circular convolution as computed by FFT.

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

• you say you get "absurd" results, what are these? – Marcus Müller Jun 9 at 20:55
• Please provide sample input and the "absurd" output that it generates – tobassist Jun 10 at 14:05
• @tobassist the sample input is included in the code. For the output, run the code. – DeltaIV Jun 10 at 14:59
• @MarcusMüller I modified the code and the question extensively, but the results still don't make sense, even now I compute the circular convolution. – DeltaIV Jun 10 at 17:24
• @tobassist I modified the question and included the output. – DeltaIV Jun 10 at 17:24

• ouch! Thanks for the pointer. Is there a numpy function to compute cyclic convolution, so that I can test the correctness of my code? – DeltaIV Jun 9 at 21:12
• implementing that in a for loop doesn't sound so bad, honestly. – Marcus Müller Jun 9 at 21:14
• do you mean implementing cyclic convolution in a for loop? – DeltaIV Jun 9 at 21:15
• A reference for circular convolution would be z = ifft(fft(x).*fft(y)) (MATLAB syntax) – Hilmar Jun 9 at 21:39