Sample compute of a DFT (Discrete Fourier Transform) with pure python code (without numpy usage)

pure_py_dft.py

#!/usr/bin/env python3
# sample compute of a DFT (Discrete Fourier Transform) with pure python code (without numpy usage)

from cmath import exp
from math import pi, sin
import random


def dft(samples):
    xl = [0.0] * len(samples)
    # for each Xl
    for l in range(len(xl)):
        # for current Xl do sum of every xs part
        for s, xs in enumerate(samples):
            xl[l] += xs * exp(-2j * pi / len(samples)) ** (l * s)
    return xl


def linspace(start, stop, n=100):
    if n < 2:
        return stop
    diff = (float(stop) - start) / (n - 1)
    return [diff * i + start for i in range(n)]


if __name__ == '__main__':
    # Number of sample points
    N = 300
    # sample spacing
    T = 0.001

    # build signal
    t = linspace(0.0, N * T, N)
    y1 = [0.0] * N

    # sin waves sum
    for i in range(N):
        # 50 Hz wave
        y1[i] += 10.0 * sin(2.0 * pi * 50 * t[i])
        # 100 Hz wave
        y1[i] += 30.0 * sin(2.0 * pi * 100 * t[i])
        # 200 Hz wave
        y1[i] += 30.0 * sin(2.0 * pi * 200 * t[i])
        # add some noise
        y1[i] += 5.0 * random.random() - 0.5

    # dft of y1
    yf1 = dft(y1)

    # format spectrum points (xf1, mf1)
    # frequency domain from 0 to 1/(2*sampling time) Hz
    xf1 = linspace(0.0, 1.0 / (2.0 * T), N // 2)
    # magnitude (absolute value of complex dft result)
    mf1 = []
    for i in range(N//2):
        mf1.append(2.0 / N * abs(yf1[i]))

    # view results with matplotlib
    import matplotlib.pyplot as plt

    plt.subplots_adjust(hspace=0.4)
    plt.subplot(211)
    plt.plot(t, y1)
    plt.xlabel("t (s)")
    plt.ylabel("y1")
    plt.subplot(212)
    plt.plot(xf1, mf1)
    plt.xlabel("xf1 (Hz)")
    plt.ylabel("mf1")
    plt.grid()
    plt.show()