1

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.])
DeltaIV
  • 111
  • 3

2 Answers2

1

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.

Marcus Müller
  • 30,525
  • 4
  • 34
  • 58
  • 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 09 '21 at 21:12
  • implementing that in a for loop doesn't sound so bad, honestly. – Marcus Müller Jun 09 '21 at 21:14
  • do you mean implementing cyclic convolution in a for loop? – DeltaIV Jun 09 '21 at 21:15
  • yes, that's what I mean. – Marcus Müller Jun 09 '21 at 21:16
  • hmmm. For a unit test I'd rather use an existing, well-tested function from a known library, than writing code myself. Anyway, can you add an implementation of circular convolution as a for loop to your question? Or should I ask a new question? – DeltaIV Jun 09 '21 at 21:20
  • seriously, if you can't trust yourself to write the convolution, you can't trust yourself to write a correct unit test, so get over it and write it :) – Marcus Müller Jun 09 '21 at 21:30
  • 2
    A reference for circular convolution would be z = ifft(fft(x).*fft(y)) (MATLAB syntax) – Hilmar Jun 09 '21 at 21:39
  • @Hilmar that's an actually useful comment: I came to the same conclusion in meantime. Just as a curiosity, do you have a link to the definition of circular convolution for finite-length sequences? Thanks! – DeltaIV Jun 10 '21 at 05:21
  • @DeltaIV um. https://en.wikipedia.org/wiki/Convolution#Circular_discrete_convolution – Marcus Müller Jun 10 '21 at 05:22
  • there was an error in my code, and it wasn't even computing the circular convolution, actually. Now I'm comparing circular convolution computed using the DHT, with circular convolution computed using the FFT. They still don't match. – DeltaIV Jun 10 '21 at 17:26
0

You forgot to roll the array after flipping. What you want is x0, x(n-1), x(n-1)... x2, x1 but you are using x(n-1),x(n-2),... x1, x0 instead.

Change

Xflip = np.flip(X)
Yflip = np.flip(Y)

to this

Xflip = np.roll(np.flip(X),1)
Yflip = np.roll(np.flip(Y),1)

to get the correct convolution.

syockit
  • 111
  • 3