1
$\begingroup$

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.])
$\endgroup$
5
  • $\begingroup$ you say you get "absurd" results, what are these? $\endgroup$ Jun 9 at 20:55
  • $\begingroup$ Please provide sample input and the "absurd" output that it generates $\endgroup$
    – tobassist
    Jun 10 at 14:05
  • $\begingroup$ @tobassist the sample input is included in the code. For the output, run the code. $\endgroup$
    – DeltaIV
    Jun 10 at 14:59
  • $\begingroup$ @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. $\endgroup$
    – DeltaIV
    Jun 10 at 17:24
  • $\begingroup$ @tobassist I modified the question and included the output. $\endgroup$
    – DeltaIV
    Jun 10 at 17:24
1
$\begingroup$

I didn't know the convolution theorem for the DHT before, but it's pretty clear that if it exists, it must be about circular convolution, just like for the DFT.

You're comparing that with acyclic convolution, so the results differing is no surprise.

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.