Inversion Labs Step functions (or square waves) are common in signal generation, switching circuits, or for timing reference. An ideal mathematical step function is a square wave that changes between the high and low state instantaneously. Although this discontinuity is impossible to achieve in physical systems it can be approximated using a Fourier series, defined as a finite sum of continuous trigonometric functions (sin and cosine).
The code below generates a continuous step function using a Fourier approximation. As an example, we'll analyze the "rise time" between the high and low (or "open" and "closed") states. Rise time is the time taken for the signal to change between specified limits, commonly 10% and 90% of output step height, and can be important for precision electronics.
This example generates a signal that switches between high (1.0) and low (-1.0) states every 5 seconds. Let's assume our application requires a rise time of less than 10 ms (0.01 seconds). The goal is to increase the number of iterations of the Fourier series until we achieve a signal that is within the rise time threshold.
First, here is a short segment of the close-to-ideal square wave:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import math
def f(t):
T = 10
if 0 < t < (T/2.):
return 1
elif t == (T/2.):
return 0
elif (T/2.) < t < T:
return -1
x = np.arange(0, 10, 0.0001) # Define the range and step-size for the output
plt.figure(figsize=(6 * 1.618, 6))
plt.plot(x, np.array([f(t) for t in x]), label='f(t)')
plt.ylim((-1.1, 1.1))
plt.xlabel('Time (s)', fontsize=12)
plt.xticks(fontsize = 12)
plt.yticks(fontsize=12)
Now, the continuous Fourier approximation, where we experiment to see how many iterations are necessary to achieve the 0.01 sec transition:
def S(t,N):
T = 10
i = 1
sum = 0.
while i <= N: # compute the Fourier series for N iterations
a = (1 / (2.*i-1.)) * math.sin((2. * (2.*i-1.) * math.pi * t) / T)
sum = sum + a
i += 1
return (4. / math.pi) * sum
x = np.arange(0, 30, 0.005) # Define the range and step-size for the output
plt.figure(figsize=(6 * 1.618, 6))
plt.plot(x, np.array([S(t, 10) for t in x]), label='S(t, 10)')
plt.plot(x, np.array([S(t, 100) for t in x]), label='S(t, 100)')
plt.plot(x, np.array([S(t, 1000) for t in x]), label='S(t, 1000)')
plt.xlabel('Time (s)', fontsize=12)
plt.ylim(-1.4, 1.9)
plt.xticks(fontsize = 12)
plt.yticks(fontsize=12)
plt.legend()
plt.show()
Note that the oscillations occurring near the transition points are called the "Gibbs Phenomenon". This is a cause of ringing artifacts in signal processing, and can be reduced with low-pass filters.
Taking a closer look at the 5-sec transition point, and plotting the ideal step function (black line), we can see that the approximation with 1000 iterations (red line) is the most accurate, with the transition between 10% and 90% of the step height occurring within 0.01 sec. Of course, within the limits of CPU capacity these approximations could be even further improved.
b = np.arange(0, 10, 0.0001)
plt.figure(figsize=(6 * 1.618, 6))
plt.plot(b, np.array([f(t) for t in b]), label='f(t)', color='k')
plt.plot(x, np.array([S(t, 10) for t in x]), label='S(t, 10)')
plt.plot(x, np.array([S(t, 100) for t in x]), label='S(t, 100)')
plt.plot(x, np.array([S(t, 1000) for t in x]), label='S(t, 1000)')
plt.axhline(y=0.8, linestyle='--', linewidth=1, color = 'k', label='10%, 90% levels')
plt.axhline(y=-0.8, linestyle='--', linewidth=1, color = 'k')
plt.xlabel('Time (s)', fontsize=12)
plt.xticks(fontsize = 12)
plt.yticks(fontsize=12)
plt.xlim(4.98, 5.02)
plt.ylim(-1.2, 1.2)
plt.legend()
plt.show()
Overall, this is a convenient piece of code to easily compute Fourier approximations, with a large variety of possible applications.
Thanks to Jesse Johnson in the Computer Science Department at University of Montana for introduction to the concepts presented in this post.