Show the code

import numpy as np
import matplotlib.pyplot as plt
from sympy import *
from IPython.display import Markdown
from sympy.solvers.ode.systems import dsolve_system
from sympy import *

0.1 Fourier Series

We begin by thinking about a string that is fixed at both ends. When a sinusoidal wave is reflected from the ends, for some frequencies the superposition of the two waves will form a standing wave with a node at each end. We characterize which standing wave is set up on the string by an integer \(n = 1, 2, 3, 4, ...\). In general n can be any positive integer, so there are in principle an infinite number of possible standing waves. (Figure 1) shows the first four possible standing waves.

Show the code

x = np.linspace(0, 1, 500) * np.pi

y1 = np.sin(x)
y2 = np.sin(2.0 * x)
y3 = np.sin(3.0 * x)
y4 = np.sin(4.0 * x)


plt.text(x[-1] + 0.1, y1[-1] + 6.4, r'$n$')


plt.plot(x, y1 + 4.4, 'r-')
plt.plot(x, -y1 + 4.4, 'b-')
plt.text(x[0] - 0.7, y1[0] + 4, 'Fundamental \n 1st Harmonic')
plt.text(x[-1] + 0.1, y1[-1] + 4.4, '1')

plt.plot(x, y2 + 2.2, 'r-')
plt.plot(x, -y2 + 2.2, 'b-')
plt.text(x[0] - 0.7, y2[0] + 2, 'First Overtone \n 2nd Harmonic')
plt.text(x[-1] + 0.1, y2[-1] + 2.2, '2')


plt.plot(x, y3 + 0.0, 'r-')
plt.plot(x, -y3 + 0.0, 'b-')
plt.text(x[0] - 0.7, y3[0] + 0, 'Second Overtone \n 3rd Harmonic')
plt.text(x[-1] + 0.1, y3[-1] + 0.2, '3')

plt.plot(x, y4 - 2.2, 'r-')
plt.plot(x, -y4 - 2.2, 'b-')
plt.text(x[0] - 0.7, y4[0] - 2.4, 'Third Overtone \n 4th Harmonic')
plt.text(x[-1] + 0.1, y4[-1] - 2.2, '4')

plt.plot([x[0], x[0]], [-4, 6], 'k-')
plt.plot([x[-1], x[-1]], [-4, 6], 'k-')

plt.axis('off')

plt.show()

For a particular standing wave, any point on the string is executing simple harmonic motion, but the amplitude depends on the position on the string. We can write the string’s displacement as a function of position and time, D(x, t), as the product of the amplitude and the harmonic oscillation. For the n-th standing wave: \[ D_n(x ,t ) =A_n (x )cos(\omega_n+\phi_n ) \qquad(1)\] If the length of the string is \(L\), then the amplitude to the \(n\)-th standing wave is: \[ A_n(x)=a_{\max}\sin(2\pi \frac{x}{\lambda_n}) \qquad(2)\] where \(\lambda_n\) is the wavelength of the standing wave, which is: \[ \lambda_n = \frac{2L}{n} \] So, for \(n = 1\) the wavelength as twice the length of the string. For \(n = 2\) the wavelength equals the length of the string, and so on. This can be easily seen in Figure 1. For a real string, such as on a guitar or violin, the amplitude \(A_\text{general}(x, t)\) will in general be pretty complicated. However, as Fourier realized long ago, that complicated vibration is just the sum of the possible standing waves. In terms of the amplitude: \[ A_\text{general}(x) = \sum_{n = 1}^\infty b_n\sin\left(2\pi \frac{x}{\lambda_n}\right) \qquad(3)\] The boundary condition that the string is fixed on both ends means that the amplitude varies as a sine function. If the string were not fixed at the ends, the boundary condition would be that at the ends of the strings the vibration will be an anti-node. In this case, the amplitude varies as a cosine function. So a more general form of Equation 3 for any boundary condition is: \[ A_\text{general}(x) = \sum_{n=1}^\infty a_n\cos\left(2\pi\frac{x}{\lambda_n}\right) + \sum_{n=1}^\infty b_n\sin\left(2\pi\frac{x}{\lambda_n}\right). \qquad(4)\]

There is one more subtlety. We have assumed that the oscillations are around \(y = 0\). A more general form would allow for oscillations about some non-zero value. Then we have: \[ A_\text{general}(x) = \frac{a_0}{2}\sum_{n=1}^\infty a_n\cos\left(2\pi\frac{x}{\lambda_n}\right) + \sum_{n=1}^\infty b_n\sin\left(2\pi\frac{x}{\lambda_n}\right). \qquad(5)\] Finally, we know that: \[ e^{-i\theta} = \cos\theta -i\sin\theta. \] so we can re-write the Fourier series in terms of complex coefficients cn as: \[ A_\text{general}(x) = \sum_{i=1}^\infty c_n e^{-i\left(2\pi\frac{x}{\lambda_n}\right)} \] For the \(n\)-th standing wave each point on the string is executing simple harmonic motion with angular frequency \(\omega_n\). Since for any wave the frequency \(f\) times the wavelength \(\lambda\) is equal to the speed of propagation of the wave down the string, we can relate the frequency and period \(T\) of the oscillation to the wavelength for the \(n\)-th standing wave: \[ \lambda_nf_n = v \] \[ f_n = \frac{\omega_n}{2\pi} = \frac{v}{\lambda_n} \] \[ \omega_n = \frac{nv}{2L} \] \[ T_n = \frac{1}{f_n} = \frac{2L}{nv} \] Each standing wave will generate a sound wave of the same frequency traveling through the air. We can write the amplitude of the sound wave as a function of time, \(y_n(t)\) , at some position away from the string as: \[ \begin{eqnarray*} y_n(t) &=& a_\max \sin\left(2\pi\frac{t}{T_n}+\phi\right) \\ &=& a_\max \sin(\omega_n t + \phi) \end{eqnarray*} \] Of course, in general the sound wave will be some complicated function of the time, but that complicated function is a sum of the sound waves from the individual standing waves on the string. Thus we can de-compose the sound wave as a function of time the same way we did for the amplitude of the standing wave as a function of position. \[ A_\text{general}(x) = \frac{a_0}{2}\sum_{n=1}^\infty a_n \cos \omega_n t + \sum_{n=1}^\infty b_n \sin \omega_n t \] Musicians call the \(n = 1\) term the fundamental: it represents the note that the string is playing. The terms with \(n > 1\) are called the overtones, and the relative amounts and phases of the overtones determines the timbre of the sound, so that a violin and a Les Paul guitar playing the same note sound different. Integrating sine and cosine functions for different values of the frequency shows that the terms in the Fourier series are orthogonal. In terms of the Kronecker delta: \[ \delta_{nm} \equiv \left\lbrace\begin{aligned}[l] 1 & m = n \\ 0 & m \neq n \end{aligned}\right. \qquad(6)\]

the orthogonality conditions are: \[ \begin{eqnarray*} \int_0^{2\pi/\omega_1} \sin (n\omega_1 t) \sin (m \omega_1 t) dt &=& \delta_{nm} \frac{\pi}{\omega_1} \\ \int_0^{2\pi/\omega_1} \cos (n\omega_1 t) \cos (m \omega_1 t) dt &=& \delta_{nm} \frac{\pi}{\omega_1} \\ \int_0^{2\pi/\omega_1} \sin (n\omega_1 t) \cos (m \omega_1 t) dt &=& 0; ~~~ \forall~~ m, n \end{eqnarray*} \qquad(7)\]

You will want to notice that the limits of the integrations in (Equation 7) are from \(0\) to \(T_1\), the period of the \(n = 1\) standing wave.

Given the function \(y_\text{general}(t)\), we can find the coefficients of the Fourier series in (Equation 6) using the orthogonality conditions to get: \[ \begin{eqnarray*} a_n &=& \frac{\omega_1}{\pi}\int_0^{2\pi/\omega_1} y_\text{general}(t)\cos(n\omega_1 t)dt \\ b_n &=& \frac{\omega_1}{\pi}\int_0^{2\pi/\omega_1} y_\text{general}(t)\sin(n\omega_1 t)dt \end{eqnarray*} \]

0.2 Fourier Transforms

In (Equation 7), the integration is over \(0\) to \(T_1\), or \(-\frac{T_1}{2}\) to \(\frac{T_1}{2}\), what about a non-periodic range? \[ Y(\omega) = \int_{-\infty}^{+\infty}y(t) e^{-i \omega t} dt \qquad(8)\] (Equation 8) defines the Fourier transform. Physically we have resolved a single pulse or wave packet \(y(t)\) into it frequency components. Notice that \(Y\) is only a function of the angular frequency, so we have transformed a function of time into a function of angular frequency. We can also define the inverse Fourier transform which takes a function of angular frequency and produces a function of time: \[ y(t) = \frac{1}{2\pi}\int_{-\infty}^{+\infty} Y(\omega)e^{+i\omega t} d\omega \qquad(9)\] If you look in other sources, you may see other conventions for the Fourier transform and the inverse Fourier transform. Some normalize the integral of Equation 8 by multiplying the integral by \(1/\sqrt{2\pi}\) and multiplying the integral in Equation 9 by the same factor of \(1/\sqrt{2\pi}\). Still other sources have the Fourier transform involve a positive exponential, with the inverse transform using the negative exponential. We will always use the conventions of Equation 8 and Equation 9 in this document. Our first example will be \(10\) oscillations of a sine wave with angular frequency \(\omega = 2.5~\text{s}^{-1}\) \[ y(t) = \left\lbrace\begin{aligned}[l] \sin 2.5 t & -4\pi \leq t \leq +4\pi \\ 0 & |t| > 4\pi \end{aligned}\right. \]

The Fourier transform is: \[ \begin{eqnarray*} Y(\omega) &=& \int_{-\infty}^{+\infty}y(t)e^{-i\omega t}dt \\ &=& \int_{-4\pi}^{+4\pi} \sin(2.5 t) e^{-i\omega t}dt \end{eqnarray*} \qquad(10)\] Since \(y(t)\) is a sine function, we expect the Fourier transform (Equation 10) to be purely imaginary.

Show the code

import scipy
from sympy import *
def y(t):
    indexes = (np.abs(t) > 4.0 * np.pi)
    result = np.sin(2.5 * t)
    result[indexes] = 0.0
    return result

it = np.linspace(-60, 60, 1200)

plt.figure(figsize=(13, 4))
plt.subplot(121)
plt.plot(it, y(it))
plt.xlabel(r'$t(s)$')
plt.ylabel(r'$y(t)$')

omega = np.linspace(-10, 10, 2000)
var('o t')
ff = integrate(sin(2.5 * t) * exp(-1.0j * o * t), (t, -4.0 * np.pi, 4.0 * np.pi))


yomega = np.zeros(omega.shape)
fff = lambdify((o), ff)
for i in range(len(omega)):
    yomega[i] = fff(omega[i]).imag
plt.subplot(122)
plt.plot(omega, yomega)
plt.xlabel(r'$\omega(t^{-1})$')
plt.ylabel(r'$Y(\omega)$')
plt.show()

Markdown(f'$${latex(ff)}$$')

First, the Fourier transform has a negative peak at \(2.5 s^{-1}\) and a positive peak at \(–2.5 s^{-1}\). The negative peak at \(+2.5 s^{-1}\) is minus the sine component of the frequency spectrum. It is negative because we chose the negative exponential for the Fourier transform, the imaginary part is minus the sine component. Signs are always a bit of a pain in Fourier transforms, especially since, as already mentioned, different implementations use different conventions.

The positive peak at \(–2.5 s^{-1}\) arises because the angular frequency could just as well have a negative value as positive one, but the sine function is anti symmetric, \(\sin\theta = -\sin(-\theta)\).

The second thing that you should notice is that there is significant overlap between the curve of the positive angular frequency solution and the negative angular frequency one.

Nonetheless, to generate a finite sine wave pulse requires the superposition of a number of frequencies in addition to the frequency of the sine wave itself.

There is also a computational issue of which you should be aware. Although it is possible to evaluate Equation 14 by hand giving a purely imaginary solution, rounding errors mean that doing it with software such as Mathematica will produce small but non-zero real terms. For this case the largest value of the calculated real component of the Fourier transform as evaluated by Mathematica is a negligible \(-5 \times 10^{-17}\). This property of software evaluation of Fourier transforms will occur again in this document.

We will now take the Fourier transform of the same \(\sin(2.5t)\) function, but this time for \(30\) oscillations.

\[ y(t) = \left\lbrace\begin{aligned}[l] \sin 2.5 t & -12\pi \leq t \leq +12\pi \\ 0 & |t| > 12\pi \end{aligned}\right. \]

The Fourier transform is: \[ \begin{eqnarray*} Y(\omega) &=& \int_{-\infty}^{+\infty}y(t)e^{-i\omega t}dt \\ &=& \int_{-12\pi}^{+12\pi} \sin(2.5 t) e^{-i\omega t}dt \end{eqnarray*} \qquad(11)\]

Show the code

import scipy
from sympy import *
def y(t):
    indexes = (np.abs(t) > 12.0 * np.pi)
    result = np.sin(2.5 * t)
    result[indexes] = 0.0
    return result

it = np.linspace(-60, 60, 1200)

plt.figure(figsize=(13, 4))
plt.subplot(121)
plt.plot(it, y(it))
plt.xlabel(r'$t(s)$')
plt.ylabel(r'$y(t)$')

omega = np.linspace(-10, 10, 2000)
var('o t')
ff = integrate(sin(2.5 * t) * exp(-1.0j * o * t), (t, -12.0 * np.pi, 12.0 * np.pi))


yomega = np.zeros(omega.shape)
fff = lambdify((o), ff)
for i in range(len(omega)):
    yomega[i] = fff(omega[i]).imag
plt.subplot(122)
plt.plot(omega, yomega)
plt.xlabel(r'$\omega(t^{-1})$')
plt.ylabel(r'$Y(\omega)$')
plt.show()

Markdown(f'$${latex(ff)}$$')

Comparing with Figure 2, you can see that the overall shape of the Fourier transform is the same, with the same peaks at \(–2.5 s^{-1}\) and \(+2.5 s^{-1}\), but the distribution is narrower, so the two peaks have less overlap. If we imagine increasing the time for which the sine wave is non-zero to the range \(-\infty\) to \(+\infty\) the width of the peaks will become zero. Physically this means that there is only one frequency component of an infinitely long sine wave pulse, and it is equal to the frequency of the sine wave itself.

In Physics, being able to resolve a signal into its frequency components is immensely useful. However, there is more Physics contained in the Fourier transform.

First, you may have already recognized the shape of the Fourier transforms in Figures 2(b) and 3(b). They are identical to the wave amplitudes of single-slit diffraction. This is not a coincidence. Although we have been thinking of the variable t as time, imagine for a moment that it is a position. Then Equation 13 is describing diffraction through a slit whose width is \(8\pi\), while Equation 15 is for a slit whose width is \(24\pi\). So the fact that the width of the distribution in Figure 3(b) is narrower than the distribution in Figure 2(b) is telling us that the width of the diffraction pattern for a narrow slit is greater than the width for a wide slit.

Here is some more Physics hidden in the Fourier transform. If we have \(N\) oscillations of a sine wave pulse of the form \(\sin(\omega t)\) , it is not too difficult to show that the width of the central maximum of the Fourier transform is: \[ \Delta \omega = \frac{\omega_0}{N} \] If we think about the photon aspect of an electromagnetic sine wave, the energy of the photon is: \[ E_\text{photon} = hf_0 = \hbar \omega_0 \] But since the frequency of the Fourier transform has angular frequencies in addition to \(\omega_0\), there is an uncertainty in the true value of the angular frequency. It is fairly reasonable to take value of the uncertainty from the width of the central maximum, Equation 17. So there is an uncertainty in the energy of the photon: \[ \begin{eqnarray*} E_\text{photon} &=& \hbar\Delta\omega \\ &=& \hbar \frac{\omega_0}{N} \end{eqnarray*} \] If we are at some position in space, the time t it takes for the wave to pass us is \(NT_0 = N 2\pi / \omega_0\). This is the uncertainty in the time when the photon actually passes us, so: \[ \Delta t_\text{photon} = 2\pi \frac{N}{\omega_0} \] The product of these two uncertainties, our uncertainty in the energy of the photon and our uncertainty about when it had that energy, is: \[ \Delta E_\text{photon} \Delta t_\text{photon} = \hbar \frac{\omega_0}{N} \times 2\pi \frac{N}{\omega_0} \] Thus the product of the uncertainties is Planck’s constant, independent of the value of the number of oscillations N or the frequency \(\omega_0\) of the sine wave. Heisenberg’s uncertainty principle actually states: \[ \Delta E_\text{photon} \Delta t_\text{photon} \geq \hbar \] We conclude that the uncertainty principle is related to the Fourier transform.

0.3 Time Series

Although theorists often deal with continuous functions, real experimental data is almost always a series of discrete data points. For 3 oscillations of the \(\sin(2.5 t)\) wave we were considering in the previous section, then, actual data might look like the dots in Figure 4. Of course, good data would include errors in at least the dependent variable if not both variables, but we will ignore that in this document. Figure 4 If we have n data points, the data will look like:

Show the code

t1 = np.linspace(0, 8, 50)
plt.plot(t1, np.sin(2.5 * t1), 'x')
t1 = np.linspace(0, 8, 500)
plt.plot(t1, np.sin(2.5 * t1), '-')
plt.show()

\[ y_0, y_1, y_2, \dots, y_{n-1} \qquad(12)\] Such data are called a time series. In a sort-of poor convention, the sampling interval is usually given the symbol \(\Delta\). For the data of Figure 4, \(\Delta =0.20 s\) and \(n = 38\). For any sampling interval, the times corresponding to the data points of Equation 12 are: \[ 0, \Delta, 2\Delta, \dots, (n-1)\Delta \qquad(13)\] For time series, we replace the integral in the Fourier transform, Equation 11, with a sum and the differential time dt with the sampling interval \(\Delta\): \[ Y_j = Y(\omega_j) = \left(\sum_{k=0}^{n-1}y_k e^{-i\omega_j t_k}\right)\times\Delta \qquad(14)\] Equation 14 is the discrete Fourier transform. As we shall soon see, we can set the value of \(\Delta\) to \(1\) without loss of generality.

The inverse discrete Fourier transform gets back the values of \(y_k\) with: \[ y_k = \frac{1}{2\pi}\left(\sum_{j=0}^{n-1}Y_je^{+i\omega_jt_k}\right)\times\delta\omega \qquad(15)\] Soon we will discuss the \(\delta\omega\) factor in Equation 15, which replaces \(d\omega\) in the continuous inverse transform Equation 10.

Counting from \(0\), as in Equation 12,Equation 13,Equation 14,Equation 15, may be new to you. In our everyday life we commonly count from \(1\), as in “one, two, three, …, \(n\)”. This is not how counting is done for time series, or many modern computer languages such as C++, Java, and Python, which count as “zero, one, two, … , \(n – 1\)”. This can take a bit of time to get used to. This counting issue is why the current year, 2011, is in the 21st century, not the 20th.

To actually use Equation 14, Equation 15, we need to know the times \(t_k\) and angular frequencies \(\omega_j\). The times are just: \[ t_k = k \Delta \qquad(16)\] The determine \(\omega_j\), we will think about Fourier series. The value of \(\omega_0 = 0\) and corresponds to a DC component in the signal. If the signal is periodic with a period \(T = n\Delta\), then the next value of the angular frequency is: \[ \omega_1 = \frac{2\pi}{T} \] The next term is: \[ \omega_2 = 2 \times \frac{2\pi}{T} \] In general: \[ \omega_j = j\frac{2\pi}{T} = j\frac{2\pi}{n\Delta} \qquad(17)\] the discrete Fourier transform Equation 14 becomes: \[ Y_j = \left(\sum_{k=0}^{n=1}y_k e^{-i2\pi\frac{jk}{n}}\right)\times \Delta \qquad(18)\] In the definition of the inverse discrete Fourier transform, Equation 15, the sum is multiplied by \(\delta\omega\), which is how much the angular frequency \(\omega_j\) changes as \(j\) goes to \(j + 1\). We have just seen that this is: \[ \delta\omega = \frac{2\pi}{T} = \frac{2\pi}{n\Delta} \qquad(19)\] So the inverse discrete Fourier transform, Equation 15, becomes: \[ y_k = \frac{1}{n}\left(\sum_{j=0}^{n-1}Y_je^{+i2\pi\frac{jk}{n}}\right)\frac{1}{\Delta} \qquad(20)\] Now, when we actually use the discrete Fourier transform, we end up with a series of values \(Y_0, Y_1, Y_2, \dots, y_{n-1}\). If we want to know the frequency of the \(Y_j\) term, we can just use Equation 17. So the factor of \(\Delta\) that multiplies the sum of Equation 18 is not needed, and we just set its value to \(1\). Similarly, the inverse discrete Fourier transform returns a series of values \(y_0, y_1, y_2, \dots, y_{n-1}\) and if we want to the know the time of the value of \(y_k\), we can just use Equation 16. So for the inverse discrete Fourier transform we can similarly just set \(\Delta = 1\). So the final form of the discrete Fourier transform is: \[ Y_j = \sum_{k=0}^{n-1}y_k e^{-i2\pi\frac{jk}{n}} \qquad(21)\] and the inverse discrete Fourier transform is: \[ y_k = \frac{1}{n}\sum_{j=0}^{n-1}Y_je^{+i2\pi\frac{jk}{n}} \qquad(22)\]

Show the code

tmax = 3.0 * 2.0 * np.pi / 2.5
t = np.arange(0, tmax, 0.2)
it = np.linspace(0, tmax, 1000)
plt.figure(figsize = (14, 4))
plt.subplot(121)
plt.plot(it, np.sin(2.5 * it), '-')
plt.plot(t, np.sin(2.5 * t), 'o')
plt.xlabel(r'$t$')
plt.ylabel(r'$y_k$')

def myfft(tarr, yarr):
    nn = len(tarr)
    res = np.zeros(tarr.shape)
    karr = np.linspace(0, nn-1, nn)
    for i in range(nn):
        res[i] = np.sum(yarr * np.exp(-2j * np.pi * i * karr / nn)).imag / nn
    return res

plt.subplot(122)
res = myfft(t, np.sin(2.5 * t))
delta = (t[2] - t[1])
idelta = 1.0 / delta
omegaj = np.zeros(t.shape)
for i in range(len(omegaj)):
    omegaj[i] = 2.0 * np.pi * i / (len(t) * delta)
plt.plot(omegaj, res * idelta, '-o')
#plt.xlim([0, 5])
plt.xlabel(r'$\omega_j$')
plt.ylabel(r'$Y_j$')
plt.show()

There may be a major surprise for you in Figure 4. You can see the negative peak, which for the continuous Fourier transforms Figure 2(b) and Figure 3(b) corresponded to the angular frequency of \(2.5 s^{-1}\). But the positive peak, corresponding to the angular frequency of \(12 –2.5s^{-1}\), is now to the far right. What has happened is that the discrete Fourier transform just returns a series of \(n = 38\) values: \[ Y_0, Y_1, Y_2, \dots, Y_{n-3}, Y_{n-2}, Y_{n-1} \] The first 19 of these correspond to the positive angular frequency values of Figure 2(b) and Figure 3(b). The other 19 values, corresponding to the negative angular frequencies, have just been appended to the end of the first 19 values.

Although the discrete Fourier transform shown in Figure 4 was evaluated with Mathematica, this way of handling the negative frequency solutions is standard for most implementations, including Python’s implementation. One reason why is that usually we are not interested in the negative frequency “alias” solution, so can just ignore it or even just throw out the last half of the Fourier transform data.

There may be another small surprise in Figure 4. The amplitude of the sampled sine wave is just 1, but the absolute value of minimum and maximum values of the transform is approximately 19, which is one-half the number of samples \(n = 38\). This is a consequence of the normalization conditions of Equation 21 and Equation 22. A symmetric normalization would multiply the sum by \(2/n\) for the discrete Fourier transform, and replace the \(1/n\) factor for the inverse transform by the same factor of \(2/n\). Using this convention, the maximum and minimum values would be what you might expect. This was also an issue for the continuous Fourier transforms of the previous section, but we ignored that in the discussion.

The negative peak in the imaginary part of the Fourier transform shown in Figure 4 occurs at the 4th value, which is \(j = 3\). From Equation 17, this corresponds to an angular frequency of: \[ \omega_3 = 3 \frac{2\pi}{n\Delta} = 3 \frac{2\pi}{3\times(0.20s)} = 2.48s^{-1} \]

Show the code

print(omegaj[3])

2.4802047265182576

It is reasonable to assume that the error in this value of one-half of the change in the value of the angular frequency from \(j\) to \(j +1\), \(\delta\omega/2\). From Equation 19 this is: \[ \frac{\delta\omega}{2} = \frac{1}{2} \frac{2\pi}{n\Delta} = 3 \frac{\pi}{3 \times (0.02s)} = 0.41s^{-1} \] So the frequency corresponding to the peak is: \[ \omega_3 = (2.48 \pm 0.41)s^{-1} \] which is well within errors of the actual frequency of the signal, \(2.5 s^{-1}\). Increasing the number of samples will reduce the uncertainty in the calculated value of the frequency.

We saw in our discussion of the continuous Fourier transform of a sine function that evaluating the integrals in software gave very small extraneous values for the real parts of the transform. This is also an issue for the discrete Fourier transform, but is compounded by the fact that the signal is sampled and not continuous. For the example of Figure 4, the maximum magnitude of the real part of the Fourier transform is about \(8\%\) of the maximum magnitude of the imaginary part. The effect of the sampling time \(\Delta\) on Fourier transforms is an immense topic, which we will not discuss here beyond saying that in general the smaller the sampling time the better.

The actual way that the discrete Fourier transforms, Equation 21, are implemented means that there is a caveat associated with the statement that the smaller the sampling time the better. If we evaluate Equation 21 using “brute force” we can define: \[ W = e^{-i\frac{2\pi}{n}} \qquad(23)\] Then the discrete Fourier transform becomes: \[ Y_j = \sum_{k=0}^{n-1}y_k e^{-i2\pi\frac{jk}{n}} = \sum_{k=0}^{n-1}y_kW^{jk} \qquad(24)\] We can think of \(y_k\) as a vector of length \(n\), and \(W\) as a matrix of dimension \(n \times k\). The multiplication of the two requires \(n^2\) calculations, and evaluating the sum requires a smaller number of operations to generate the powers of \(W\). Thus the number of calculations necessary to evaluate the Fourier transform is proportional to \(n^2\). Doubling the number of points in the time series quadruples the time necessary to calculate the transform, and tripling the number of points requires nine times as much time. For large data sets, then, the time necessary to calculate the discrete Fourier transform can become very large.

However, there is a brilliant alternative way of doing the calculation that is was reinvented by Cooley and Tukey in 1965. It is called the fast Fourier transform. The idea is that we split the sum into two parts: \[ Y_j = \sum_{k=0}^{n/2-1}y_{2k}e^{-i2\pi\frac{j2k}{n}} + \sum_{k=0}^{n/2-1}y_{2k+1}e^{-i2\pi} \qquad(25)\] The first sum involves the even terms \(y_{2k}\), and the second one the odd terms \(y_{2k + 1}\). Using the definition of \(W\) in Equation 23 we can write this as: \[ \begin{eqnarray*} Y_j &=& Y_j^{k~even} + W^k\sum_{k=0}^{n/2-1}y_{2k+1}e^{-i2\pi\frac{jk}{n/2}} \\ &=&Y_j^{k~even} + W^k Y_j^{k~odd} \end{eqnarray*} \qquad(26)\] We can apply the same procedure recursively. Eventually, if \(n\) is a power of two, we end up with no summations at all, just a product of terms.

An analogy to the algorithm of the fast Fourier transform is a method to determine the number of hairs on your head. Just counting all the hairs would be a very long process. However, imagine that you divide your scalp into two equal pieces. Then the total number of hairs on your head is \(2\) times the number of hairs in one of the two pieces. If you divide one of those pieces in half, the total number of hairs on your head is \(2^2\) times the number in that sample. Dividing that piece in half means that the total is \(2^3\) times the number in the new sample. If you keep dividing the area of your scalp in half a total of M times, then eventually you get down to a small enough piece that you can easily count the number of hairs in it, and the total number of hairs is \(2^M\) times the hairs in the sample area. In the limit where you divide the areas of your scalp enough times that the sample contains just one hair, then the number of hairs on your head is just \(2^M\).

It turns out that the number of calculations required using the fast Fourier transform is proportional to \(n\log_2n\). So doubling the number of points in the time series only doubles the time necessary to do the calculation, and tripling the number of points increases the time by about \(4.75\). This is a big win over the brute force method of doing the calculation.

Any competent implementation of the fast Fourier transform does not require that the number of data points in the times series be a power of two, but if not it will need to use some brute force calculations at least at the end. In one test, a time series of \(e^{-t/100}\) was generated for \(t\) from \(0\) to \(1000.001\) with \(\Delta=0.001\). The number of data points was \(n = 1 000 001\), and in one computing environment Mathematica took \(0.89 s\) to calculate the Fourier transform. The value of the last data point is \(e^{1000/100} = 0.0000454\), which is nearly zero. \(48 575\) zeroes were appended to the dataset, so the total length became \(1 048 576 = 20^2\). Mathematica took 0.20 s to calculate the Fourier transform of this larger dataset, which is over four times faster. 15

Although speed of calculation is not an issue for the small dataset of \(38\) samples of a sine wave that we are considering here, the lesson to be learned is that for large datasets for which you wish to use the fast Fourier transform, you should design the experiment so that the number of samples is a power of \(2\).

0.4 Python’s Implementation

The Python programming language has an implementation of the fast Fourier transform in its scipy and numpy library. Below we will write a single program, but will introduce it a few lines at a time.

Show the code

tmin = 0
tmax = 2.4 * np.pi
delta = 0.2 
t = np.arange(tmin, tmax, delta)
y = np.sin(2.5 * t)
plt.figure(figsize=(14, 4))
plt.subplot(121)
plt.plot(t, y, 'bo')
plt.title('y(t) = sin(2.5 t)')
plt.xlabel('t (s)')
plt.ylabel('y(t)')

Y = np.fft.fft(y)
ff= np.fft.fftfreq(len(y), t[2]-t[1])
omega = 2.0 * np.pi * ff
plt.subplot(122)
plt.plot(omega, Y.imag)
plt.show()

0.5 FFT solver for Poission equation

Show the code

import feyn

from sklearn.datasets import load_diabetes
import pandas as pd

from feyn.tools import split

# Load diabetes dataset into a pandas dataframe
dataset = load_diabetes()
df_diabetes = pd.DataFrame(dataset.data, columns=dataset.feature_names)
df_diabetes['response'] = dataset.target

# Train/test split
train, test = split(df_diabetes, ratio=[0.6, 0.4], random_state=42)

# Instantiate a QLattice
ql = feyn.QLattice(random_seed=42)

models = ql.auto_run(
    data=train,
    output_name='response'
)
# Select the best Model
best = models[0]
best

Show the code

# Plot simplified model with only 'bp' manually fixed
best.plot_response_1d(
    data = train,
    by = 'bmi',
    input_constraints = {'bp': train.bmi.quantile(0.25)}
)

1 1. Generating the signal

We will be using a cosine wave of 3 Hz frequency, as the input signal.

Show the code

import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns
sns.set()

Show the code

f = 3.0
t = np.arange(0,2,0.001) 

cos_wave = np.cos(2*np.pi*f*t)
# cos_wave = 2*m.cos(2*np.pi*f*t) + 5*m.cos(2*np.pi*f*2*t)

plt.rcParams["figure.figsize"] = (12,4)
plt.plot(t,cos_wave)
plt.title("Cosine Wave with 3 Hz frequency")
plt.ylabel("Amplitude")
plt.xlabel('Time (in seconds)');

2 2. Warping the input signal wave around

Just to summarise the idea put forth in the video, we need to try warping the input signal graph for a number of sampling frequencies and then by keeping the track of the x-coordinate of the center of mass, considering the warped graph as a piece made of metal wire, we can estimate the frequencies present in the input signal.
The range for frequency search is kept here as 0 to 10 Hz.

Show the code

r_cord = []
min_freq_range = 0.0
max_freq_range = 10.0
sf_list = np.arange(min_freq_range, max_freq_range, 0.1)
for sf in sf_list:
    r_cord.append( [(cos_wave[i], t[i]*sf*2*np.pi) for i in range(len(t)) ] )

Show the code

x_cord , y_cord = [], []
for l in range(len(r_cord)):
    x_cord.append( [amp*np.cos(theta) for (amp,theta) in r_cord[l]] )
    y_cord.append( [amp*np.sin(theta) for (amp,theta) in r_cord[l]] )

Below cell plots all the beautiful curves that pop up while doing this for each value of the sampling frequency in our list.
Will take a little bit of time to execute, because the number of plots may be high, and once every plot is made, only then the whole figure pops up.
Totally worth the wait !! The Red dot denotes the Center of Mass of the graph.

Show the code

mean_list = []

plt.rcParams["figure.figsize"] = (15,110)
for l in range(len(r_cord)):
    plt.subplot(int(len(r_cord)/4)+1, 4, int(l+1))
    plt.plot(x_cord[l], y_cord[l])
    plt.plot(np.mean(x_cord[l]), np.mean(y_cord[l]), 'or' )
    plt.title("Sampling Freq = "+str(round(sf_list[l], 2))+" Hz")
    
    # Storing the COM for plotting later
    x_mean = np.sum(x_cord[l])
    mean_list.append(x_mean)

3 3. Plotting COM v/s Sampling Frequency

Show the code

plt.rcParams["figure.figsize"] = (12,4)
plt.xlabel("Frequeny (in Hz)")
plt.xticks(np.arange(min(sf_list), max(sf_list), 0.5))
sns.set()
plt.plot(sf_list,mean_list);

For clear estimation of Frequencies present, we need to smoothen the graph.
Below is the smoothened version of the graph, with a clear peak near around 3 Hz.

Show the code

plt.rcParams["figure.figsize"] = (12,4)
smoothed = [i if i>0 and i>0.2*max(mean_list) else 0 for i in mean_list]
plt.plot(sf_list, smoothed)
plt.xlabel("Frequeny (in Hz)")
plt.xticks(np.arange(min(sf_list), max(sf_list), 0.5));

Plotting the same graph as a Bar chart, as the Frequencies are depicted in a typical Fourier Analysis

Show the code

plt.rcParams["figure.figsize"] = (12,4)
plt.xticks(np.arange(min(sf_list), max(sf_list), 0.5))
plt.xlabel("Frequeny (in Hz)")
plt.bar(sf_list, smoothed, width = 0.06);