Fourier transform

We begin with the usual loading of modules.

import numpy as np # for numerics
import sympy as sp # for symbolics
import matplotlib.pyplot as plt # for plots!

Let’s consider a periodic function \(f\) with period \(T\) (T). Each period, the function has a triangular pulse of width \(\delta\) (pulse_width) and height \(\delta/2\).

period = 15 # period
pulse_width = 2 # pulse width

First, we plot the function \(f\) in the time domain. Let’s begin by defining \(f\).

def pulse_train(t,T,pulse_width):
  f = lambda x:pulse_width/2-abs(x) # pulse
  tm = np.mod(t,T)
  if tm <= pulse_width/2:
    return f(tm)
  elif tm >= T-pulse_width/2:
    return f(-(tm-T))
  else:
    return 0

Now, we develop a numerical array in time to plot \(f\).

N = 151 # number of points to plot
tpp = np.linspace(-period/2,5*period/2,N) # time values
fpp = np.array(np.zeros(tpp.shape))
for i,t_now in enumerate(tpp):
  fpp[i] = pulse_train(t_now,period,pulse_width)

Now we plot.

fig, ax = plt.subplots()
ax.plot(tpp,fpp,'b-',linewidth=2) # plot
plt.xlabel('time (s)')
plt.xlim([-period/2,3*period/2])
plt.xticks(
  [0,period],
  [0,'$T='+str(period)+'$ s']
)
plt.yticks([0,pulse_width/2],['0','$\delta/2$'])
plt.draw()

Consider the following argument. Just as a Fourier series is a frequency domain representation of a periodic signal, a Fourier transform is a frequency domain representation of an aperiodic signal (we will rigorously define it in a moment). The Fourier series components will have an analog, then, in the Fourier transform. Recall that they can be computed by integrating over a period of the signal. If we increase that period infinitely, the function is effectively aperiodic. The result (within a scaling factor) will be the Fourier transform analog of the Fourier series components.

Let us approach this understanding by actually computing the Fourier series components for increasing period \(T\) using . We’ll use sympy to compute the Fourier series cosine and sine components \(a_n\) and \(b_n\) for component \(n\) (n) and period \(T\) (T).

x, a_0, a_n, b_n = sp.symbols('x, a_0, a_n, b_n', real=True)
delta, T = sp.symbols('delta, T', positive=True)
n = sp.symbols('n', nonnegative=True)
an = sp.integrate(
  2/T*(delta/2-sp.Abs(x))*sp.cos(2*sp.pi*n/T*x),
  (x,-delta/2, delta/2) # otherwise zero
).simplify()
bn = 2/T*sp.integrate(
  (delta/2-sp.Abs(x))*sp.sin(2*sp.pi*n/T*x),
  (x, -delta/2, delta/2) # otherwise zero
).simplify()
print(sp.Eq(a_n,an), sp.Eq(b_n,bn))

Eq(a_n, Piecewise((T*(1 - cos(pi*delta*n/T))/(pi**2*n**2), n > 0), (delta**2/(2*T), True))) Eq(b_n, 0)

Furthermore, let us compute the harmonic amplitude
(f_harmonic_amplitude): $$\begin{align} C_n = \sqrt{a_n^2 + b_n^2} \end{align}$$ which we have also scaled by a factor \(T/\delta\) in order to plot it with a convenient scale.

C_n = sp.symbols('C_n', positive=True)
cn = sp.sqrt(an**2+bn**2)
print(sp.Eq(C_n, cn))

Eq(C_n, Piecewise((T*Abs(cos(pi*delta*n/T) - 1)/(pi**2*n**2), n > 0), (delta**2/(2*T), True)))

Now we lambdify the symbolic expression for a numpy function.

cn_f = sp.lambdify((n, T, delta), cn)

Now we can plot. Write a function to plot pulses in the time domain with the corresponding frequency spectrum.

def plot_pulses_and_spectrum(T, pulse_width, omega_max):
  n_max = round(omega_max*T/(2*np.pi)) # max harmonic
  n_a = np.linspace(0,n_max,n_max+1)
  omega = 2*np.pi*n_a/T
  fig, ax = plt.subplots(1, 2)
  plt.sca(ax[0])
  for i in range(0, 3):
    tpp = np.linspace(-T/2, 5*T/2,N)
    fpp = np.array(np.zeros(tpp.shape))
    for i,t_now in enumerate(tpp):
      fpp[i] = pulse_train(t_now, T, pulse_width)
    plt.plot(tpp, fpp, 'b-', linewidth=2)
  plt.xlim([-T/2, 3*T/2])
  plt.xticks([0, T], [0, '$T='+str(T)+'$ s'])
  plt.yticks([0, pulse_width/2], ['0', '$\delta/2$'])
  plt.xlabel('time (s)')
  plt.sca(ax[1])
  plt.stem(
    omega, cn_f(n_a, T, pulse_width)*T/pulse_width, 'bo-'
  )
  plt.xlim([0, omega_max])
  plt.ylim([0, 1.1])
  plt.xlabel('Frequency $\omega$ (rad/s)')
  plt.ylabel('$C_n T/\delta$')
  return fig

Now we plot the pulses and their spectra for \(T \in [5,15,25]\) rad/s and \(\delta = 2\).

omega_max = 12  # Maximum frequency to plot
fig = plot_pulses_and_spectrum(5, pulse_width, omega_max)
plt.draw()

fig = plot_pulses_and_spectrum(15, pulse_width, omega_max)
plt.draw()

fig = plot_pulses_and_spectrum(25, pulse_width, omega_max)
plt.draw()

The line spectra are shown in the right-hand columns of the plots above. Note that with our chosen scaling, as \(T\) increases, the line spectra reveal a distinct waveform.

Let \(F\) be the continuous function of angular frequency \(\omega\) $$\begin{align} F(\omega) = \frac{\delta}{2} \cdot \frac{\sin^2(\omega \delta/4)}{(\omega \delta/4)^2}. \end{align}$$ First, we plot it.

def F(w):
  return pulse_width/2*np.sin(w*pulse_width/4)**2 / \
    (w*pulse_width/4)**2
N = 201 # number of points to plot
wpp = np.linspace(0.0001, omega_max,N)
Fpp = []
for i in range(0,N):
    Fpp.append(F(wpp[i])) # build array of function values
fig, ax = plt.subplots()
plt.plot(wpp, Fpp, 'b-', linewidth=2) # plot
plt.xlim([0, omega_max])
plt.yticks([0, pulse_width/2],['0','$\delta/2$'])
plt.xlabel('Frequency $\omega$ (rad/s)')
plt.ylabel('$F(\omega)$')
plt.show()

The plot of \(F\) is clearly emerging from the preceding line spectra as the period \(T\) increases.

Now we are ready to define the Fourier transform and its inverse. We will define the Fourier transform in two ways: as a trigonometric transform and as a complex transform. We begin with the trigonometric transform and its inverse.

The Fourier transform is a generalization of the Fourier series to aperiodic functions (i.e., functions with infinite period). The complex form of the Fourier transform is more convenient for analysis and computation, as we will see.

So now we have defined the Fourier transform. There are many applications, including solving differential equations and frequency domain
representations—called spectra—of time domain functions.

There is a striking similarity between the Fourier transform and the Laplace transform, with which you are already acquainted. In fact, the Fourier transform is a special case of a Laplace transform with Laplace transform variable \(s = j\omega\) instead of having some real component. Both transforms convert differential equations to algebraic equations, which can be solved and inversely transformed to find time-domain solutions. The Laplace transform is especially important to use when an input function to a differential equation is not absolutely integrable and the Fourier transform is undefined (for example, our definition will yield a transform for neither the unit step nor the unit ramp functions). However, the Laplace transform is also preferred for initial value problems due to its convenient way of handling them. The two transforms are equally useful for solving steady state problems. Although the Laplace transform has many advantages, for spectral considerations, the Fourier transform is the only game in town.

A table of Fourier transforms and their properties can be found in section B.2.

Example 7.3

Consider the aperiodic signal y(t) = u_s(t)e^−at with u_s the unit step function and a > 0. The signal is plotted below. Derive the complex frequency spectrum and plot its magnitude and phase.

The signal is aperiodic, so the Fourier transform can be computed from [@eq:fourier_transform]: $$\begin{align*} Y(\omega) &= \int_{-\infty}^\infty y(t) e^{j\omega t} dt \\ &= \int_{-\infty}^\infty u_s(t) e^{-a t} e^{j\omega t} dt \tag{def.\ of $y$} \\ &= \int_{0}^\infty e^{-a t} e^{j\omega t} dt \tag{$u_s$ effect} \\ &= \int_{0}^\infty e^{(-a + j\omega) t} dt \tag{multiply} \\ &= \left.\frac{1}{-a+j\omega} e^{(-a + j\omega) t}\right|_{0}^\infty dt \tag{antiderivative} \\ &= \frac{1}{-a+j\omega} \left( \lim_{t\rightarrow \infty} e^{(-a + j\omega) t} - e^0 \right) \tag{evaluate} \\ &= \frac{1}{-a+j\omega} \left( \lim_{t\rightarrow \infty} e^{-a t} e^{j\omega t} - 1 \right) \tag{arrange} \\ &= \frac{1}{-a+j\omega} \left( (0) (\text{complex with mag} \le 1) - 1 \right) \tag{limit} \\ &= \frac{-1}{-a+j\omega} \tag{consequence} \\ &= \frac{1}{a-j\omega} \\ &= \frac{a+j\omega}{a+j\omega}\cdot\frac{1}{a-j\omega} \tag{rationalize} \\ &= \frac{a+j\omega}{a^2 + \omega^2}. \end{align*}$$

The magnitude and phase of this complex function are straightforward to compute: $$\begin{align*} |Y(\omega)| &= \sqrt{\Re(Y(\omega))^2 + \Im(Y(\omega))^2} \\ &= \frac{1}{a^2 + \omega^2}\sqrt{a^2 + \omega^2} \\ &= \frac{1}{\sqrt{a^2 + \omega^2}} \\ \angle Y(\omega) &= \arctan(\omega/a). \end{align*}$$

Now we can plot these functions of ω. Setting a = 1 (arbitrarily), we obtain the plots of figure 7.11.