3 Fourier Series and Fourier Transform
Learning Objectives
After reading this chapter, you should be able to:
3.1 Motivation
In the previous chapter we looked at the response of a single degree of freedom system excited by a sinusoidal force. We know that eventually when we want to compute the motions of a ship or offshore structure in the ocean, the excitation due to waves and wind will not be sinusoidal but will be random. Thus, there is a need to understand how we can scale from the simple case of sinusoidal forcing to a random forcing. In this chapter we will move towards this by learning about Fourier series and Fourier transforms.
3.2 Fourier series
In the previous chapter we looked at a sinusoidal forcing of the type \(F(t) = F_0 cos(\omega t)\). This forcing as seen in Figure 3.1 will repeat itself over time and has a period of \(T = 2\pi/\omega\).
Now let’s consider a funtion \(F(t)\) that is continuous and periodic over the time interval \([-\frac{T}{2}, \frac{T}{2}]\). An example is shown in Figure 3.2, where the function is periodic with a period \(T=2\pi\) but is not sinusoidal.
Fourier showed that as long as the function \(F(t)\) is smooth, continuous and periodic, it can be expressed as a trignometric series:
\[\begin{align} F(t) = \frac{a_0}{2} + \sum_{k=1}^{\infty} a_k \cos(k\omega_0 t) + b_k \sin(k\omega_0 t) \label{eq-real-fourier-series} \end{align}\]
where \(\omega_0 = 2\pi/T\) is known as the fundametal frequency of the Fourier series and \(T\) is the fundametal period of the function.
We will now compare this exercise of expressing a function as a Fourier series with the more familiar concept of projecting a vector in 3D vector space along its basis vectors. We know that in 3D space, a vector \(\vec{\boldsymbol{r}}\) is expressed in terms of its projection along the basis vectors as:
\[\begin{align} \vec{\boldsymbol{r}} = a \boldsymbol{\hat{i}} + b \boldsymbol{\hat{j}} + c \boldsymbol{\hat{k}} \end{align}\]
where \(\boldsymbol{\hat{i}}\), \(\boldsymbol{\hat{j}}\) and \(\boldsymbol{\hat{k}}\) are unit vectors along the cardinal directions (\(x\), \(y\) and \(z\) directions) and the values \(a\), \(b\) and \(c\) are the projection of the vector \(\vec{\boldsymbol{r}}\) along \(\boldsymbol{\hat{i}}\), \(\boldsymbol{\hat{j}}\) and \(\boldsymbol{\hat{k}}\).
We can think of a Fourier series as projecting the funtion \(F(t)\) along the basis functions \(\cos(k\omega_0 t)\) and \(\sin(k\omega_0 t)\) for \(k=1,2,...,\infty\). The only difference here is that we are projecting functions into function space and unlike the 3D vector space that has only 3 basis vectors, here we have infinte number of basis functions. The coefficients \(a_k\) and \(b_k\) denote the projection of the function \(F(t)\) along \(\cos(k\omega_0 t)\) and \(\sin(k\omega_0 t)\) respectively.
In the vector space, the cardinal vectors \(\boldsymbol{\hat{i}}\), \(\boldsymbol{\hat{j}}\) and \(\boldsymbol{\hat{k}}\) are orthogonal to each other. This means that when we project one cardinal vector along another, the projection is zero. This idea of projecting one vector along another is denoted by the dot product (or scalar product) operation and the orthogonality is expressed mathematically as:
\[\begin{align} \boldsymbol{\hat{i}}.\boldsymbol{\hat{j}} = \boldsymbol{\hat{i}}.\boldsymbol{\hat{k}} = \boldsymbol{\hat{j}}.\boldsymbol{\hat{k}} = 0 \label{eq-vec-orthogonality} \end{align}\]
The idea of dot product also transfers over to the function spaces. The dot product between two functions \(f(t)\) and \(g(t)\) with domain \([-\frac{T}{2}, \frac{T}{2}]\) is defined as:
\[\begin{align} <f,g> = \int_{-T/2}^{T/2} f(t) g^{*}(t) dt \end{align}\]
where \(g^{*}(t)\) is the complex conjugate of \(g(t)\). If the functions are real valued, then the dot product reduces to:
\[\begin{align} <f,g> = \int_{-T/2}^{T/2} f(t) g(t) dt \end{align}\]
Just as we had the orthogonality of basis vectors as shown in \(\eqref{eq-vec-orthogonality}\), the basis functions \(\cos(k\omega_0 t)\) and \(\sin(k\omega_0 t)\) for \(k=1,2,...,\infty\) are also orthogonal to each other.
\[\begin{align} <\cos(k\omega_0 t), \cos(m\omega_0 t)> &= \int_{-T/2}^{T/2} \cos(k\omega_0 t) \cos(m\omega_0 t) dt \nonumber \\ &= \int_{-T/2}^{T/2} \frac{1}{2}\left\{\cos((k+m)\omega_0 t) + \cos((k-m)\omega_0 t)\right\} dt \nonumber \\ &= \begin{cases} 0 & \text{for } k \neq m \\ \frac{T}{2} & \text{for } k = m \end{cases} \end{align}\]
\[\begin{align} <\sin(k\omega_0 t), \sin(m\omega_0 t)> &= \int_{-T/2}^{T/2} \sin(k\omega_0 t) \sin(m\omega_0 t) dt \nonumber \\ &= \int_{-T/2}^{T/2} \frac{1}{2}\left\{\cos((k-m)\omega_0 t) - \cos((k+m)\omega_0 t)\right\} dt \nonumber \\ &= \begin{cases} 0 & \text{for } k \neq m \\ \frac{T}{2} & \text{for } k = m \end{cases} \end{align}\]
\[\begin{align} <\cos(k\omega_0 t), \sin(m\omega_0 t)> &= \int_{-T/2}^{T/2} \cos(k\omega_0 t) \sin(m\omega_0 t) dt \nonumber \\ &= \int_{-T/2}^{T/2} \frac{1}{2}\left\{\sin((k+m)\omega_0 t) - \sin((k-m)\omega_0 t)\right\} dt \nonumber \\ &= 0 \end{align}\]
Using the orthogonality, the coefficients of the Fourier series can be found as:
\[\begin{align} a_0 &= \frac{1}{T} <F, \cos(0 \times \omega_0 t)> = \frac{1}{T} <F, 1> \\ a_k &= \frac{2}{T} <F, \cos(k \omega_0 t)> \\ b_k &= \frac{2}{T} <F, \sin(k \omega_0 t)> \end{align}\]
Figure 3.3 shows the Fourier series approximation of the non-sinusoidal function shown in Figure 3.2 as we progressively increase the number of terms considered for summation. It can be seen that with 10 terms the match is fairly good.
3.3 Complex Fourier series
For a real valued function \(F(t)\) defiend over \([-\frac{T}{2}, \frac{T}{2}]\), the Fourier series is given by \(\eqref{eq-real-fourier-series}\). This can further be expressed in terms of complex basis functions as:
\[\begin{align} F(t) &= \frac{a_0}{2} + \sum_{k=1}^{\infty} a_k \cos(k\omega_0 t) + b_k \sin(k\omega_0 t) \nonumber \\ &= \frac{a_0}{2} + \sum_{k=1}^{\infty} \frac{1}{2} \left[(a_k - ib_k)e^{ik\omega_0 t} + (a_k + ib_k)e^{-ik\omega_0 t} \right] \nonumber \\ &= \sum_{k=-\infty}^{\infty} c_k e^{ik\omega_0 t} \end{align}\]
where:
\[\begin{align} c_k = \begin{cases} \frac{1}{2} a_0 & k = 0 \\ \frac{1}{2} (a_k - ib_k) & k > 0 \\ \frac{1}{2} (a_k + ib_k) & k < 0 \end{cases} \end{align}\]
The same principle of orthogonality also holds for the complex basis functions too as shown below:
\[\begin{align} <e^{ik\omega_0 t}, e^{im\omega_0 t}> &= \int_{-T/2}^{T/2} e^{ik\omega_0 t} e^{-im\omega_0 t} dt \nonumber \\ &= \begin{cases} 0 & k \neq m \\ T & k = m \end{cases} \end{align}\]
Using the orthogonality, the coefficients of the Fourier series can be found as:
\[\begin{align} c_k = \frac{1}{T}<F, e^{ik\omega_0 t}> \end{align}\]
It is interesting to note that for a real function \(F(t)\), the complex Fourier coefficients have a symmetry as shown below:
\[\begin{align} c_{-k} = c_k^{*} \end{align}\]
3.4 Fourier transform
The Fourier series we have discussed so far for functions periodic over a finite domain \([-\frac{T}{2}, \frac{T}{2}]\). Note that the consecutive term frequencies are discretely separated by the fundamental frequency \(\omega_0\). The Fourier series can be written as:
\[\begin{align} F(t) &= \sum_{k=-\infty}^{\infty} c_k e^{ik\omega_0 t} \nonumber \\ &= \sum_{k=-\infty}^{\infty} \frac{1}{T}<F, e^{ik\omega_0 t}> e^{ik\omega_0 t} \nonumber \\ &= \sum_{k=-\infty}^{\infty} \frac{\omega_0}{2\pi} \int_{-T/2}^{T/2} F(\xi)e^{-ik\omega_0 \xi} d\xi ~e^{ik\omega_0 t} \end{align}\]
We now extend this idea of Fourier series to functions defined over the domain \((-\infty, \infty)\) by taking the limit \(T \rightarrow \infty\). Note that this means that the fundamental frequency \(\omega_0 \rightarrow d\omega\) and the sum over a discrete set of frequencies (each consecutive one separated by \(\omega_0\)) can be considered as an integral over continuous set of frequencies \(\omega\) in the limit.
\[\begin{align} F(t) &= \lim_{\substack{T\rightarrow \infty \\ \omega_0 \rightarrow d\omega}} \sum_{k=-\infty}^{\infty} \frac{\omega_0}{2\pi} \int_{-T/2}^{T/2} F(\xi)e^{-ik\omega_0 \xi} d\xi ~e^{ik\omega_0 t} \nonumber \\ &= \frac{1}{2\pi} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} F(\xi)e^{-i\omega \xi} d\xi ~e^{i\omega t} d\omega \end{align}\]
The inner integral is the Fourier transform of the continuous function \(F(t)\) defined over the domain \((-\infty, \infty)\).
\[\begin{align} \hat{F}(\omega) = \mathcal{F}(F(t)) = \int_{-\infty}^{\infty} F(t)e^{-i\omega t} dt \end{align}\]
Similarly the outer integral is the inverse Fourier transform of the continuous function \(\hat{F}(\omega)\) defined over the domain \((-\infty, \infty)\).
\[\begin{align} F(t) = \mathcal{F}^{-1}(\hat{F}(\omega)) = \frac{1}{2\pi} \int_{-\infty}^{\infty} \hat{F}(\omega) e^{i\omega t} d\omega \end{align}\]
The two integrals are known as Fourier transform pair. Both integrals converge as long as \(\int_{-\infty}^{\infty} |F(t)| dt < \infty\) and \(\int_{-\infty}^{\infty} |\hat{F}(\omega)| d\omega < \infty\). The Fourier transform and inverse Fourier transform are both linear operators and accommodate the principle of superposition.
\[\begin{align} \mathcal{F}(af(t) + bg(t)) &= a \mathcal{F}(f) + b \mathcal{F}(g) = a \hat{f}(\omega) + b \hat{f}(\omega) \\ \mathcal{F}^{-1}(a\hat{f}(\omega) + b\hat{g}(\omega)) &= a \mathcal{F}^{-1}(\hat{f}) + b \mathcal{F}^{-1}(\hat{g}) = a f(t) + b f(t) \end{align}\]
The Fourier transform gives us another way to compute the derivative of a function.
\[\begin{align} \frac{dF}{dt} &= \frac{d}{dt} \left(\frac{1}{2\pi}\int_{-\infty}^{\infty} \hat{F}(\omega)e^{i\omega t} d\omega\right) \nonumber\\ &= \frac{1}{2\pi}\int_{-\infty}^{\infty} \hat{F}(\omega)\frac{d}{dt}\left(e^{i\omega t}\right) d\omega \nonumber\\ &= \frac{1}{2\pi}\int_{-\infty}^{\infty} \hat{F}(\omega) i\omega e^{i\omega t} d\omega \nonumber\\ &= \mathcal{F}^{-1}(i\omega \hat{F}(\omega)) \end{align}\]
Thus the Fourier transform of a derivative is equivalent to mulitplying the Fourier transform of the function with \(i\omega\). This usefult property allows converting partial differntial equations into ordinary differential equations. We will use this property later when we discuss hydrodynamics of floating structures.
Another useful property of the Fourier transform is that it converts the convolution operation between two functions in the time domain into a product operation in the frequency domain. The convolution operation between two functions \(f(t)\) and \(g(t)\) is defined as:
\[\begin{align} f * g (t) = \int_{-\infty}^{\infty} f(\tau)g(t - \tau) d\tau \end{align}\]
The Fourier transform of the convolution operation yields:
\[\begin{align} \mathcal{F}(f * g) &= \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f(\tau)g(t - \tau) d\tau e^{-i\omega t} dt \nonumber \\ &= \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f(\tau)g(t - \tau) e^{-i\omega (t-\tau)} e^{-i\omega \tau}d\tau dt \nonumber \\ &= \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f(\tau)g(t - \tau) e^{-i\omega (t-\tau)} dt e^{-i\omega \tau}d\tau \nonumber \\ &= \int_{-\infty}^{\infty} f(\tau) \int_{-\infty}^{\infty} g(\xi) e^{-i\omega \xi} d\xi e^{-i\omega \tau}d\tau \nonumber \\ &= \int_{-\infty}^{\infty} f(\tau) \hat{g}(\omega) e^{-i\omega \tau}d\tau \nonumber \\ &= \int_{-\infty}^{\infty} f(\tau) e^{-i\omega \tau}d\tau \hat{g}(\omega)\nonumber \\ &= \hat{f}(\omega) \hat{g}(\omega) \end{align}\]
Table 3.1 lists some commonly used Fourier transform pairs. Here \(u(t)\) is the unit step function (\(u(t) = 1\) for \(t \geq 0\) and \(0\) otherwise), \(\text{sgn}(t)\) is the signum function (\(\text{sgn}(t) = 1\) for \(t > 0\), \(-1\) for \(t < 0\)) and \(\delta(t)\) is the Dirac delta function.
| Sl. No. | \(f(t)\) | \(\hat{f}(\omega)\) |
|---|---|---|
| 1 | \(\delta(t)\) | \(1\) |
| 2 | \(1\) | \(2\pi\delta(\omega)\) |
| 3 | \(\delta(t - t_0)\) | \(e^{-i\omega t_0}\) |
| 4 | \(e^{i\omega_0 t}\) | \(2\pi\delta(\omega - \omega_0)\) |
| 5 | \(\cos(\omega_0 t)\) | \(\pi\left[\delta(\omega - \omega_0) + \delta(\omega + \omega_0)\right]\) |
| 6 | \(\sin(\omega_0 t)\) | \(\dfrac{\pi}{i}\left[\delta(\omega - \omega_0) - \delta(\omega + \omega_0)\right]\) |
| 7 | \(u(t)\) | \(\pi\delta(\omega) + \dfrac{1}{i\omega}\) |
| 8 | \(\text{sgn}(t)\) | \(\dfrac{2}{i\omega}\) |
| 9 | \(e^{-\alpha t}u(t)\), \(\alpha > 0\) | \(\dfrac{1}{\alpha + i\omega}\) |
| 10 | \(t\,e^{-\alpha t}u(t)\), \(\alpha > 0\) | \(\dfrac{1}{(\alpha + i\omega)^2}\) |
| 11 | \(e^{-\alpha \lvert t \rvert}\), \(\alpha > 0\) | \(\dfrac{2\alpha}{\alpha^2 + \omega^2}\) |
| 12 | \(u(t)\cos(\omega_0 t)\) | \(\dfrac{\pi}{2}\left[\delta(\omega - \omega_0) + \delta(\omega + \omega_0)\right] + \dfrac{i\omega}{\omega_0^2 - \omega^2}\) |
| 13 | \(u(t)\sin(\omega_0 t)\) | \(\dfrac{\pi}{2i}\left[\delta(\omega - \omega_0) - \delta(\omega + \omega_0)\right] + \dfrac{\omega_0}{\omega_0^2 - \omega^2}\) |
| 14 | \(u(t)e^{-\alpha t}\cos(\omega_0 t)\), \(\alpha > 0\) | \(\dfrac{\alpha + i\omega}{(\alpha + i\omega)^2 + \omega_0^2}\) |
| 15 | \(u(t)e^{-\alpha t}\sin(\omega_0 t)\), \(\alpha > 0\) | \(\dfrac{\omega_0}{(\alpha + i\omega)^2 + \omega_0^2}\) |
| 16 | \(f(t) = 1\) for \(\lvert t \rvert < \tau/2\), \(0\) otherwise | \(\dfrac{2\sin(\omega\tau/2)}{\omega}\) |
| 17 | \(f(t) = \dfrac{\sin(Wt)}{\pi t}\) | \(\hat{f}(\omega) = 1\) for \(\lvert \omega \rvert < W\), \(0\) otherwise |
| 18 | \(f(t) = \tau - \lvert t \rvert\) for \(\lvert t \rvert < \tau\), \(0\) otherwise | \(\dfrac{4\sin^2(\omega\tau/2)}{\omega^2}\) |
| 19 | \(e^{-\alpha t^2}\), \(\alpha > 0\) | \(\sqrt{\dfrac{\pi}{\alpha}}\, e^{-\omega^2/4\alpha}\) |
| 20 | \(\displaystyle\sum_{n=-\infty}^{\infty} c_n e^{in\omega_0 t}\) | \(\displaystyle 2\pi\sum_{n=-\infty}^{\infty} c_n \delta(\omega - n\omega_0)\) |
The last entry shows that the Fourier transform of a periodic signal, expressed as a complex Fourier series (as in the previous section), is a train of impulses in the frequency domain located at the harmonics \(n\omega_0\) and weighted by the corresponding Fourier series coefficients \(c_n\).
3.5 Application to single degree of freedom system
Let’s revisit our differential equation in the forced vibration case:
\[\begin{align} \frac{d^2x}{dt^2} + 2\zeta\omega_n\frac{dx}{dt} + \omega_n^2 x = f(t) = \frac{F(t)}{m} \end{align}\]
Taking a Fourier transform of our differential equation and using the properties of Fourier transform for the derivatives yields:
\[\begin{align} \mathcal{F}(f(t)) &= \mathcal{F} \left(\frac{d^2x}{dt^2} + 2\zeta\omega_n\frac{dx}{dt} + \omega_n^2 x\right) \nonumber \\ \hat{f}(\omega) &= \left((i\omega)^2 + 2\zeta\omega_n (i\omega) + \omega_n^2 \right) \hat{x}(\omega) \nonumber\\ \frac{\hat{x}(\omega)}{\hat{f}(\omega) / \omega_n^2} &= \frac{1}{\left(1 - \frac{\omega^2}{\omega_n^2}\right) + i2\zeta\frac{\omega}{\omega_n}} = H(\omega) \equiv \hat{h} (\omega) \end{align}\]
the same RAO \(\hat{h} (\omega) \equiv H(\omega)\) that we saw in the previous chapter. Thus the Fourier transform of the response \(\hat{x}(\omega)\) is given by:
\[\begin{align} \hat{x}(\omega) = \frac{\hat{h}(\omega) \hat{f}(\omega)}{\omega_n^2} = \frac{\hat{h}(\omega) \hat{F}(\omega)}{k} \end{align}\]
Now we can apply the inverse Fourier transform and use the property of the Fourier transform for the convolution integral to get:
\[\begin{align} x(t) &= \mathcal{F}^{-1}(\hat{x}(\omega)) = \mathcal{F}^{-1}\left(\frac{\hat{h}(\omega) \hat{F}(\omega)}{k}\right) \nonumber \\ &= \frac{1}{k} \int_{-\infty}^{t} F(t) h(t - \tau) d\tau \label{eq-irf-conv} \end{align}\]
where \(h(t) = \mathcal{F}^{-1}(\hat{h}(\omega))\) is known as the impulse response function (IRF) for the system and is the inverse Fourier transform of the RAO. Note the integral upper limit is upto \(t\) and not \(\infty\). This is so to make sure causality is maintained. In physical systems, a force (cause) at a specific time can affect the response (effect) only at that time and in the future. We cannot have a force (cause) in the future affect the response (effect) in the past. Such systems where this cause and effect relationship is maintained are called causal systems. If the upper limit of the integral in \(\eqref{eq-irf-conv}\) is allowed to be \(\infty\), then a force in the future \(F(t')\) with \(t' > t\) can affect the response \(x(t)\) at time \(t\). This breaks causality and is unphysical.
Using the Fourier transform pairs shown in Table 3.1, we can compute the closed form expression for the impulse response function as:
\[\begin{align} h(t) &= \mathcal{F}^{-1}\left(\frac{1}{\left(1 - \frac{\omega^2}{\omega_n^2}\right) + i2\zeta\frac{\omega}{\omega_n}}\right) \nonumber \\ &= \mathcal{F}^{-1}\left(\frac{\omega_n^2}{\left(\omega_n^2 - \omega^2\right) + i2\zeta\omega_n\omega}\right) \nonumber \\ &= \mathcal{F}^{-1}\left(\frac{\omega_n^2}{\left(\zeta\omega_n + i\omega\right)^2 + \omega_n^2(1-\zeta^2)}\right) \nonumber \\ &= \mathcal{F}^{-1}\left(\frac{\omega_n^2}{\omega_n\sqrt{1 - \zeta^2}}\frac{\omega_d}{\left(\zeta\omega_n + i\omega\right)^2 + \omega_d^2}\right) \nonumber \\ &= \frac{\omega_n}{\sqrt{1 - \zeta^2}} e^{-\zeta\omega_nt}\sin(\omega_d t) \end{align}\]
The plot of the impulse response function for \(\omega_n = 3.0\) rad/s and \(\zeta=0.1\) is shown in Figure 3.4.
Note that this response shown in Figure 3.4 is the response of the mass when the input excitation is \(F(t) = k\delta(t)\) where \(\delta(t)\) is the Dirac delta function or the impulse function. The Dirac delta is zero everywhere except where the argument of the function is zero (in this case \(t=0\)).
Notion of sum of impluses
A general force time history can be thought of as an infinite series of impulses of strength \(F(\tau)\) at each time \(\tau\). The total response \(x(t)\) at some time \(t\) is then the sum of the individual responses to each of these impulses summed at time \(t\) due to all the impulses that came prior to time \(t\). Response of one impulse that occured at time \(\tau\) is \(F(\tau)h(t-\tau)/k\) and affects the total response for all \(t>\tau\). We need to sum across all impulses that happended at all times \(\tau < t\). This is given by the following (which turns out to be the convolution integral):
\[\begin{align} x(t) = \frac{1}{k} \int_{-\infty}^{t} F(\tau)h(t-\tau) d\tau \label{eq-convolution-irf} \end{align}\]
Notice that this a huge improvement as compared to our previous chapter. There we were restricted to only sinusoidal forcing. But now for any irregular forcing \(F(t)\) we can compute the response \(x(t)\) of our single degree of freedom system.
In general computing the convolution term shown in \(\eqref{eq-convolution-irf}\) is computationally expensive as we need to slide one signal across the other and for each sliding offset we multiply the functions and calculate the area under it. The same process in Fourier domain is a lot cheaper as it is simply a multiplication of the Fourier transfrom of the functions. In adopting the computation in the Fourier domain we are essentially separating the irregular force time history into its component regular waves (through Fourier transform), computing the response of each of these regular wave components and summing them together (inverse Fourier transform) to get the response to the irregular force.