Synthetic signals

import numpy as np
import matplotlib.pyplot as plt
import scipy

def f0(x):
    return np.mod(np.floor(2*x), 2)

# Alternate version of f0
def f0a(x):
    return ((x - np.floor(x)) >= 0.5) * 1
    
# Alternate version of f0
def f0b(x):
    return (1 - (-1) ** np.floor(2*x)) / 2

def f1(x):
    return x - np.floor(x)

def f2(x):
    return np.abs(x - np.round(x))

# Alternate version of f2
def f2a(x):
    return np.minimum(x - np.floor(x),np.ceil(x) - x)

def f3(x, a = 0.25):
    return 1 * (x - np.floor(x) < a)

def f4(x):
    return np.abs(np.sin(np.pi*x))

functions = [f0, f1, f2, f3, f4]
function_names = ["Square wave", "Sawtooth wave", "Triangle wave", "Pulse wave", "Rectified sine wave"]

ts = np.linspace(0.0001, 3.999, 1013)
fig, ax = plt.subplots(5, 1, figsize=(10, 10))
fig.tight_layout(pad=3.0)
for i in range(5):
    ax[i].plot(ts, functions[i](ts))
    ax[i].set_title(function_names[i])

for i in range(5):
    plt.figure(figsize=(10, 2))
    plt.plot(ts, functions[i](ts))
    plt.title(function_names[i])
    plt.show()

The cell below checks that the functions f0a and f0b are really the same as f0, and that f2a is really the same as f2.

print(np.max(np.abs(f0(ts) - f0a(ts))))
print(np.max(np.abs(f0(ts) - f0b(ts))))
print(np.max(np.abs(f2(ts) - f2a(ts))))

0.0
0.0
0.0

The array p0 is the same as the discrete Fourier transform of the vector f0(ts) divided by N. The array q0 contains the first N values of the continuous Fourier transform of the function f0. When k is small relative to N, the values p0[k] and q0[k] are close. The corresponding facts are also true for p1, p2, p3 and p4.

N = 1024
ns = np.arange(N)
ns_pos = np.arange(1,N)
ns_odd = np.arange(1,N,2)
ns_even = np.arange(0,N,2)

theta = 2 * np.pi / N
omega = np.exp(-1j * theta)

p0 = np.zeros(N, dtype=complex)
p0[0] = 1/2
p0[1::2] = 2 / N / (omega ** ns_odd - 1)

q0 = np.zeros(N, dtype=complex)
q0[0] = 1/2
q0[1::2] = 1j / np.pi / ns_odd

p1 = np.zeros(N, dtype=complex)
p1[0] = (1 - 1/N)/2
p1[1:] = 1 / N / (omega ** ns_pos - 1)

q1 = np.zeros(N, dtype=complex)
q1[0] = 1/2
q1[1:] = 1j / (2 * np.pi) / ns_pos

p2 = np.zeros(N, dtype=complex)
p2[0] = 1/4
p2[1::2] = 2 / (N * N) / (np.cos(ns_odd * theta) - 1)

q2 = np.zeros(N, dtype=complex)
q2[0] = 1/4
q2[1::2] = -1 / (np.pi ** 2) / (ns_odd ** 2)

p3 = np.zeros(N, dtype=complex)
p3[0] = 1/4
# p3[1:] = (1 - omega ** (256 * np.arange(1,1024))) / (1 - omega ** (np.arange(1,1024)))
p3[1:] = (1 - (-1j) ** ns_pos) / (1 - omega ** ns_pos) / N

q3 = np.zeros(N, dtype=complex)
q3[0] = 1/4
q3[1:] = -1j / (2 * np.pi) * (1 - (-1j) ** ns_pos) / ns_pos

p4 = np.zeros(N, dtype=complex)
p4[0] = 1 / N / np.tan(theta/4)
p4[1:] = np.sin(theta/2) / (np.cos(np.arange(1,N)*theta) - np.cos(theta/2)) / N

q4 = 2 / np.pi / (1 - 4 * ns ** 2)

ps = [p0, p1, p2, p3, p4]
qs = [q0, q1, q2, q3, q4]

We can check the above claims as follows. The first tuple printed below is $(0, 11.746, 1.952)$. The $0$ means that the line refers to f0. The $11.746$ means that the array p0 agrees with the discrete Fourier transform as computed by scipy.fft.fft to nearly $12$ digits of accuracy (the base 10 logarithm of the error is $-11.746$). The $1.952$ means that the base $10$ logarithm of the difference between the first quarter of p0 and the first quarter of q0 is $-1.952$.

def check_all():
    ts = np.arange(0, 1, 1./N)
    for i in range(5):
        p = ps[i]
        q = qs[i]
        u = scipy.fft.fft(functions[i](ts)) / N
        p_err = -np.round(np.log(np.linalg.norm(p - u)) / np.log(10),3)
        q_err = -np.round(np.log(np.linalg.norm((q - u)[:N//4])) / np.log(10),3)
        print((i, p_err, q_err))

check_all()

(0, 11.746, 1.952)
(1, 12.033, 2.102)
(2, 12.593, 5.425)
(3, 11.873, 1.952)
(4, 11.946, 5.079)

def show_function(i):
    ts = np.arange(0, 2, 1./N)
    f = functions[i]
    p = ps[i]
    q = qs[i]
    fig, ax = plt.subplots(1, 2, figsize=(15, 5))
    ax[0].plot(ts, f(ts))
    ax[0].set_title(function_names[i])
    ax[1].plot(np.abs(p)[:64], 'r.', label="p")
    ax[1].plot(np.abs(q)[:64], 'g+', label="q")
    ax[1].legend()

show_function(0)

show_function(1)

show_function(2)

show_function(3)

show_function(4)