hartley_transforms.py

import math


def fht(vector):
    n = len(vector)
    if n <= 1:
        return vector[:]
    elif n == 2:
        return [vector[0] + vector[1], vector[0] - vector[1]]

    v1 = fht(vector[::2]) * 2
    v2 = fht(vector[1::2]) * 2

    result = [0] * n

    cos = 1
    sin = 0
    cosw = math.cos(2 * math.pi / n)
    sinw = math.sin(2 * math.pi / n)
    for i in range(n):
        result[i] = v1[i] + v2[i]*cos + v2[-i]*sin

        cos, sin = cos * cosw - sin * sinw, sin * cosw + cos * sinw

    return result


def cyclic_convolution(a, b):
    n = len(a)
    a = fht(a)
    b = fht(b)

    result = [0] * n
    for i in range(n):
        # if b[ i] = c + s, then
        #    b[-i] = c - s

        e = (b[i] + b[-i])/2  # (c+s) + (c-s) = 2c  # "Even part" = cosine component
        o = (b[i] - b[-i])/2  # (c+s) - (c-s) = 2s  # "Odd part" = sine component
        result[i] = (a[i] * e + a[-i] * o) / n

    return fht(result)


def slow_cyclic_convolution(a, b):
    n = len(a)
    c = [0] * n

    for i in range(n):
        for j in range(n):
            c[(i + j) % n] += a[i] * b[j]

    return c


if __name__ == "__main__":
    n = int(input("n = "))

    a = [0] * n

    while True:
        f = input("f = ")
        if f == "":
            break
        d = float(input("a[f] += "))
        a[int(f)] += d

    a = fht(a)

    import matplotlib.pyplot as plt
    plt.plot(a)
    plt.show()

    if n <= 16:
        for x in a:
            print(round(x, 5))