Discrete Fourier Transform From Scratch

dft.py

# X_k = \sum_{n=0}^{N-1} x_n \cdot e^{-\frac{j 2 \pi k n}{N}}
# X_k = \sum_{n=0}^{N-1} x_n \cdot [\cos(-\frac{2 \pi k n}{N}) + j\sin(-\frac{2 \pi k n}{N})]

import math
import cmath

def generate_data(timespan=64):
    data = []
    for i in range(timespan):
        ni = i/timespan
        data.append(math.sin(2 * math.pi * 3 * ni + math.pi/6) + 2*math.sin(2 * math.pi * 7 * ni + math.pi/6))
    return data

def dft(signal):
    output = []
    N = len(signal)

    for k in range(N):
        acc = 0
        for n in range(N):
            acc += signal[n] * cmath.exp(-(1j * 2 * math.pi * k * n) / N)
        output.append(acc)

    return output

def idft(fourier_series):
    output = []
    N = len(fourier_series)

    for k in range(N):
        acc = 0
        for n in range(N):
            acc += fourier_series[n] * cmath.exp((1j * 2 * math.pi * k * n) / N)
        output.append((1/N) * acc)

    return output

input_data = generate_data()
# for i in input_data:
#     print(i)

fourier = dft(input_data)
freqs = [abs(x)/(len(input_data)/2) for x in fourier]
phases = [math.atan2(x.imag, x.real) for x in fourier]
# print(fft)

# inverse = idft(fourier)
# for i in inverse:
#     print(i.real)

for i in range(len(freqs)):
    f = freqs[i]
    p = phases[i]
    if abs(f) < 1e-10:
        continue
    print(f"{i}Hz: {f} with {p} rad phase")

# Known good!
# fft_res = numpy.fft.fft(input_data)
# print(fft_res)#.round(15))
# print(numpy.abs(fft_res))#.round(15))
# print(numpy.angle(fft_res))