Fourier Series and Frequency Response – Fusion of Engineering, Control, Coding, Machine Learning, and Science

In this post, we introduce frequency responses of linear systems and provide a brief introduction to Fourier series. Most importantly, we perform a real physical experiment of observing a frequency response of a resistor-capacitor (RC) circuit. A video about this post is given below. A YouTube video accompanying this post is given below. The video contains the frequency response measurement experiment.

Fourier Series

First, we provide a brief introduction to Fourier Series. Mr. Fourier, who was a French mathematician, claimed that any periodic function (even a periodic function with square corners!) can be used to be mathematically expressed as a sum of sinusoids! Sounds truly amazing.

So, let $f(t)$ be a periodic function with the period of $T$ . According to Fourier, we can represent this function as follows

(1) $\begin{align*}f(t)=\sum_{m=-\infty}^{m=+\infty} c_{m} e^{jm\omega_{0}t}\end{align*}$

where $\omega_{0}= 2\pi /T$ is the \textit{fundamental frequency}, $c_{m} \in \mathbb{C}$ is a complex coefficient, $j$ is an imaginary unit. So how to find the Fourier Series expansion \eqref{FourierSeries} for a given function $f(t)$ . It can be shown that the coefficients are defined as follows:

(2) $\begin{align*}c_{m}=\frac{1}{T} \int_{T} f(t) e^{-jm\omega_{0}t} dt\end{align*}$

Let us determine the Fouries series expansion for a periodic function shown in figure below.

The period of this function is $T=2$ . This is a square wave function ( or a periodic pulse function) that is often encountered in signal analysis and control engineering. In its period $[0,T]$ This function can be formally expressed as follows:

(3) $\begin{align*}& f(t)=A, \;\; t \in (0,T/2) \notag \\& f(t)=0,\;\; t \in (T/2,T) \notag \end{align*}$

where $A$ is the function amplitude.
So let us apply the formula \eqref{fouriesCoefficients}

(4) $\begin{align*}& c_{m}=\frac{1}{T} \int_{0}^{T} f(t) e^{-jm\omega_{0} t} dt = \frac{1}{T} \int_{0}^{T/2} 1 e^{-jm\frac{2\pi}{T} t} dt \notag \\ & c_{m}=\frac{A}{(-jm2\pi)} \big(e^{-jm\pi}-1\big) \notag \\& c_{m}=\frac{A}{m2\pi}\Big(sin(m\pi)-2\frac{j}{2}\big(1-cos(m\pi)\big) \Big)\end{align*}$

Next, we use the following trigonometric identities:

(5) $\begin{align*}sin^{2}(\frac{x}{2})=\frac{1}{2}\big(1-cos(x) \big) \notag \\sin(x)=2sin(\frac{x}{2}) cos(\frac{x}{2}) \notag \end{align*}$

By using these identities in the previous expression, we obtain

(6) $\begin{align*}c_{m}=\frac{A}{m\pi}sin(\frac{m\pi}{2})\Big(cos(\frac{m\pi}{2})-jsin(\frac{m\pi}{2}) \Big)\end{align*}$

Next, by using the Euler’s formula: $e^{xj}=cos(x)+jsin(x)$ , we obtain

(7) $\begin{align*}c_{m}= \frac{A}{2} \frac{ sin(\frac{m\pi}{2}) }{\frac{m\pi}{2}}e^{-j\frac{m\pi}{2}}\end{align*}$

Now, we have to resolve possible indeterminate cases (when $m=0$ ). We notice that for $m=0$ , $c_{m}$ takes the form of $0/0$ . We can write:

(8) $\begin{align*}c_{0}=\lim_{m->0} \frac{A}{2} \frac{ sin(\frac{m\pi}{2}) }{\frac{m\pi}{2}}e^{-j\frac{m\pi}{2}} = \lim_{m->0} \frac{A}{2} \frac{sin(\frac{m\pi}{2}) }{\frac{m\pi}{2}} = \frac{A}{2} \end{align*}$

where we used L’Hôpital’s rule to compute the limit (we can also use the limit of $sin(x)/x$ to obtain the last result). On the other hand, we have:

(9) $\begin{align*}c_{k}e^{jk\omega_{0}t}= \frac{A}{2} \frac{sin(\frac{k\pi}{2}) }{\frac{k\pi}{2}} e^{j\big(k\omega_{0}t-\frac{k\pi}{2} \big) } \notag \\ c_{-k}e^{-jk\omega_{0}t}= \frac{A}{2} \frac{sin(\frac{k\pi}{2}) }{\frac{k\pi}{2}} e^{-j\big(k\omega_{0}t-\frac{k\pi}{2} \big) } \notag \end{align*}$

By substituting the Euler’s formula, and by substituting the last expression as well as the expression for $c_{0}$ given in \eqref{c0} in the expression \eqref{ FourierSeries }, we obtain

(10) $\begin{align*}f(t)\approx f_{k_{max}}= \frac{A}{2}+\sum_{k=1}^{k_{max}} A \frac{sin(\frac{k\pi}{2})}{ \frac{k\pi}{2}} cos(k\omega_{0}t-\frac{k\pi}{2})\end{align*}$

The following figure demonstrates the performance of the Fourier series expansion.

# -*- coding: utf-8 -*-
"""

Fourier series expansion - analytical

"""
import numpy as np
import matplotlib.pyplot as plt
import math
A=1
T=2

x=np.arange(-T,2*T,0.01)
y=np.zeros(shape=(x.shape[0],))
for i in range(x.shape[0]):
    if (x[i]<(-T/2)) or ((x[i]>0) and (x[i]<T/2)) or ((x[i]>T) and (x[i]<3*T/2)):
        y[i]=1
    

fig= plt.figure(figsize=(8,4))
plt.xlim(-2.5, 4.5)
plt.ylim(-0.2, 1.2)
plt.axis('equal')
plt.ylabel('f(t)')
plt.xlabel('time - t')
plt.plot(x,y,'b', linewidth=2.5, label='Periodic function')
plt.savefig('periodic_fuction.eps')


k_max=100
w0=2*math.pi/T
yapprox=[]

for i in range(x.shape[0]):
    sum=0
    for k in np.arange(1,k_max+1):
        sum=sum+A*((math.sin(k*math.pi/2))/(k*math.pi/2))*math.cos(k*w0*x[i]-k*math.pi/2)
    yapprox.append(sum+0.5*A)   
    

fig= plt.figure(figsize=(8,4))
plt.xlim(-2.5, 4.5)
plt.ylim(-0.2, 1.2)
plt.axis('equal')
plt.ylabel('f(t)')
plt.xlabel('time - t')
plt.plot(x,y,'b', linewidth=2.5,  label='Periodic function')
plt.plot(x,yapprox,'r', linewidth=2.5,  label='Fourier series approximation')
plt.legend()
plt.savefig('approximationk_100.eps')

Frequency Response

Consider a linear dynamical system S shown in Fig. 1.

Figure 1: A block diagram of a dynamical system.

The system input is denoted by $u(t)$ and the system output is $y(t)$ . The input is an external signal affecting the system dynamics, and the output is a signal that we usually observe and whose dynamical behavior is of interest to us.

There are several ways of mathematically describing the relationship between the inputs and outputs. We can use differential equations to relate inputs and outputs. Differential equations can be used to argue about system stability and by solving them we can compute the system output for a given input. However, in many situations, we would like to have a simpler representation of the system behavior, that can be more useful from an engineering standpoint.

Frequency response is such a system representation. First, we will motivate the need for frequency response functions.

We know that (almost) every (periodic) function can be represented as a sum of harmonic functions (sin() and cos() functions). This follows from the Fourier series expansion. For example, consider a periodic ramp function:

(11) $\begin{align*}f(x)=\frac{x}{\pi},\;\; -\pi<x <\pi\end{align*}$

and where the period is $2\pi$ . Its Fourier series expansion is:

(12) $\begin{align*}f(x)=\frac{2}{\pi} \sum_{i=0}^{\infty} \frac{(-1)^{i+1}}{i} \sin(ix)\end{align*}$

Motivated by this, let us write the input function $u(t)$ as follows:

(13) $\begin{align*}u(t)=\sum_{i=0}^{\infty} a_{i} \cos(\omega_{i} t)\end{align*}$

where $\omega_{i}$ is the angular frequency and $a_{i}$ is a real constant. Since the system is linear, the sytem response to the input $u(t)$ is

(14) $\begin{align*}y(t)=\sum_{i=0}^{\infty} a_{i} y_{i}\end{align*}$

where $y_{i}$ is the system response to the harmonic function $\cos(\omega_{i} t)$ . So, the system response can be computed by summing the system responses to basic cos() functions.

This motivates us to consider the following problem. What is the system response when the input function is $\cos(\omega t)$ ? The answer is that it is again a cos() function, with unaltered angular frequency and with phase and amplitude changes that depend on the frequency. How can we show this?

First of all, we know that the system response can be computed as a convolution of an impulse response and the input:

(15) $\begin{align*}y(t)=\int_{-\infty}^{\infty} h(\tau) u(t-\tau) d\tau\end{align*}$

where $h(t)$ is the impulse response. The impulse response is the system response for the impulse function $\delta (t)$ (Dirac delta function). This function is zero everywhere except at 0 where it has an infinite value. Furthermore, its integral from minus infinity to plus infinity is 1. From the engineering point of view, this function is obtained by exciting the system with a low duration and a high-intensity pulse. Let us assume that the input is given as follows:

(16) $\begin{align*}u(t)=e^{st}\end{align*}$

where $s$ is generally a complex number. Substituting \eqref{ inputU} in \eqref{convolution1}, we obtain:

(17) $\begin{align*}y(t) = \int_{-\infty}^{\infty} h(\tau) e^{s(t-\tau)} d\tau \notag \\ y(t) = \underbrace{\int_{-\infty}^{\infty} h(\tau) \e^{-s\tau} d\tau}_{H(s)} \cdot e^{st} \end{align*}$

where the term

(18) $\begin{align*}H(s)= \int_{-\infty}^{\infty} h(\tau) e^{-s\tau} d\tau \end{align*}$

is the system transfer function. Notice that the system transfter function is the Laplace transformation of the impulse response function. Let us now use \eqref{transfer1} to compute the system response to a cos() function. Using the Euler’s relation, we can write:

(19) $\begin{align*}u(t)=A_{u} \cos(\omega t)= \frac{A_{u}}{2}\big(e^{j\omega t}+ e^{-j\omega t} \big)=u_{1}(t)+u_{2}(t)\end{align*}$

where $j$ is the imaginary unit and

(20) $\begin{align*}u_{1}(t)=\frac{A_{u}}{2} e^{j\omega t} ,\;\; u_{2}(t)=\frac{A_{u}}{2} e^{-j\omega t} \end{align*}$

Since the system is linear, the reponse to the input $u(t)$ in \eqref{cosEuler} is a sum of the responses to $u_{1}(t)$ and $u_{2}(t)$ . Let us compute the response to $u_{1}$ . Taking into account that $s=j\omega$ and using \eqref{transfer1}, we can write:

(21) $\begin{align*}y_{1}(t)=\frac{A_{u}}{2}H(j\omega)e^{j\omega t}\end{align*}$

On the other hand, $H(j\omega)$ is a complex number and it can be represented in the polar form: $H(j\omega)= M(\omega)e^{j \phi (\omega)}$ , where $M(\omega)$ is the magnitude and $\phi (\omega)$ is the phase. Using this transformation, \eqref{y1First} can be expressed as follows:

(22) $\begin{align*}y_{1}(t)=\frac{A_{u}M(\omega)}{2}e^{j(\omega t+\phi(\omega))}\end{align*}$

For $u_{2}(t)$ , we can write:

(23) $\begin{align*}y_{2}(t)= \frac{A_{u}}{2}H(-j\omega)e^{-j\omega t} \end{align*}$

Since $H(-j\omega) =M(\omega)e^{-j\phi(\omega)}$ , we have:

(24) $\begin{align*} y_{2}(t)= \frac{A_{u}M(\omega)}{2}e^{-j(\omega t+\phi(\omega))} \end{align*}$

Finally, we have:

(25) $\begin{align*}y(t)=y_{1}(t)+y_{2}(t)=M(\omega) A_{u} \frac{1}{2} \big(e^{j(\omega t +\phi(\omega))}+ e^{-j(\omega t +\phi(\omega))} \big)\end{align*}$

From the Euler’s relation, it follows that the second term in the last expression is equal to $\cos(\omega t+\phi(\omega))$ . That is,

(26) $\begin{align*} y(t)=\underbrace{M(\omega) A_{u}}_{A_{o}} \cos(\omega t+\phi(\omega)) \end{align*}$

To conclude, if the system input is $A_{u}\cos(\omega t)$ , the output function $M(\omega) A_{u} \cos(\omega t+\phi(\omega))$ . The conclusion is that the system does not alter the frequncy, but it generally speaking alters the amplitude and the phase.

The frequency response is characterized by the following two functions:

(27) $\begin{align*} M(\omega)=\frac{A_{o}}{ A_{u} } \\\phi(\omega) \label{frequencyResponse2} \end{align*}$

The first one in \eqref{frequencyResponse1} is the magnitude response. It tells us the amplification of the amplitude of the input signal. The second one in \eqref{frequencyResponse2} is the phase response. The phase response gives us information about the delay of the output signal with respect to the input signal.

In the sequel, we derive a frequency response of an RC circuit, and we present experimental results of observing the response for several frequency values. The RC circuit is shown in Fig. 2.

A differential equation describing the relation between the input ( $V_{i}$ ) and output ( $V_{o}$ ) voltages is

(28) $\begin{align*}T \dot{V}_{o}(t)+V_{o}(t)=V_{i}(t) \end{align*}$

where $T=RC$ is a time constant. Taking the Laplace transform of the last equation we obtain:

(29) $\begin{align*}T\big(s V_{o}(s)-V_{o}(0) \big)+V_{o}(s)=V_{i}(s)\end{align*}$

where $V_{o}(s)$ and $V_{s}(s)$ are the Laplace transforms of $V_{o}(t)$ and $V_{i}(t)$ , respectively, and $V_{0}(0)$ . For the computation of the transfer function, we assume that the initial conditions are zero. Consequently, we obtain:

(30) $\begin{align*}H(s)=\frac{V_{o}(s)}{V_{i}(s)}=\frac{1}{Ts+1}\end{align*}$

For $s=j\omega$ , we have

(31) $\begin{align*}H(j\omega)= \frac{1}{jT\omega +1}=\frac{1}{\sqrt{1+T^2\omega^2}} e^{-\arctan(T\omega)}\end{align*}$

From the last equation we conclude:

(32) $\begin{align*}M(\omega)= \frac{1}{\sqrt{1+T^2\omega^2}}, \;\; \phi(\omega) = -\arctan(T\omega) \end{align*}$

Next, we simulate the frequency response. The Python code is given below:

<pre class="wp-block-syntaxhighlighter-code">import numpy as np
import matplotlib.pyplot as plt

R=550
C=100*10**(-9)
T=R*C

w=0.5*np.arange(1000000)
M=(np.sqrt(1+(T*w)**(2)))**(-1)
phi=-np.arctan(T*w)*(180/np.pi)

plt.plot(w,M)
plt.title('Magnitude response')
plt.xlabel('<img src="https://aleksandarhaber.com/wp-content/ql-cache/quicklatex.com-707fcec15e450815730425a6607a1858_l3.png" class="ql-img-inline-formula quicklatex-auto-format" alt="&#92;&#111;&#109;&#101;&#103;&#97;" title="Rendered by QuickLaTeX.com" height="8" width="11" style="vertical-align: 0px;"/>')
plt.ylabel('M(<img src="https://aleksandarhaber.com/wp-content/ql-cache/quicklatex.com-707fcec15e450815730425a6607a1858_l3.png" class="ql-img-inline-formula quicklatex-auto-format" alt="&#92;&#111;&#109;&#101;&#103;&#97;" title="Rendered by QuickLaTeX.com" height="8" width="11" style="vertical-align: 0px;"/>)')
plt.savefig('magnitude.eps')
plt.show()


plt.plot(w, phi)
plt.title('Phase response')
plt.xlabel('<img src="https://aleksandarhaber.com/wp-content/ql-cache/quicklatex.com-707fcec15e450815730425a6607a1858_l3.png" class="ql-img-inline-formula quicklatex-auto-format" alt="&#92;&#111;&#109;&#101;&#103;&#97;" title="Rendered by QuickLaTeX.com" height="8" width="11" style="vertical-align: 0px;"/>')
plt.ylabel('<img src="https://aleksandarhaber.com/wp-content/ql-cache/quicklatex.com-5b2be26c0c1341f54b29baddda771346_l3.png" class="ql-img-inline-formula quicklatex-auto-format" alt="&#92;&#112;&#104;&#105;" title="Rendered by QuickLaTeX.com" height="16" width="11" style="vertical-align: -4px;"/>(<img src="https://aleksandarhaber.com/wp-content/ql-cache/quicklatex.com-707fcec15e450815730425a6607a1858_l3.png" class="ql-img-inline-formula quicklatex-auto-format" alt="&#92;&#111;&#109;&#101;&#103;&#97;" title="Rendered by QuickLaTeX.com" height="8" width="11" style="vertical-align: 0px;"/>) [deg]')
plt.savefig('phase.eps')
plt.show()</pre>

Here are the magnitude and phase responses plots:

Another important concept we need to introduce is the concept of the cutoff frequency. The cutoff frequency is usually defined as the frequency $\omega$ at which the magnitude response is equal to $\sqrt{2}/2=1/\sqrt{2}\approx 0.707$ . That is, at this frequency the amplitude of the output signal is $0.707$ times the amplitude of the input signal. Let us see what is the cutoff frequency of the RC circuit. We have

(33) $\begin{align*}M(\omega_{c})=\frac{1}{\sqrt{1+T^{2}\omega^{2}_{c}}}=\frac{1}{\sqrt{2}}\end{align*}$

From here, we have:

(34) $\begin{align*}\omega_{c}=\frac{1}{T}=\frac{1}{RC}\end{align*}$

The experimental results of measuring certain points of an RC circuit are given in the Figure below. The results are generated for $R=330\; [\Omega]$ and $C=484 10^{-9} \; [nF]$ .

The code for generating these plots is given below and for more details see the video.

<pre class="wp-block-syntaxhighlighter-code">import numpy as np
import matplotlib.pyplot as plt
import math
R=330
C=484*10**(-9)
T=R*C
freq_theoretical=1*np.arange(1500)
M=1/(np.sqrt(1+(T*(2*np.pi*freq_theoretical))**(2)))
phi=-np.arctan(T*(2*np.pi*freq_theoretical))*(180/np.pi)

cutoff=1/(R*C*2*math.pi)

frequency = np.array([195, 399, 590,  804, 992, 1196])
vout_max   = np.array([4.56, 4.24, 3.8, 3.36, 3.04,  2.72])
vin_max  = np.array([4.56, 4.48,   4.40, 4.24, 4.24,  4.16])
M_measured= vout_max/vin_max


plt.plot(freq_theoretical,M,label='theoretical')
plt.plot(frequency,M_measured,'rx',markersize=10,label='experimental')
plt.title('Magnitude response')
plt.xlabel('frequency [Hz]')
plt.ylabel('M(<img src="https://aleksandarhaber.com/wp-content/ql-cache/quicklatex.com-707fcec15e450815730425a6607a1858_l3.png" class="ql-img-inline-formula quicklatex-auto-format" alt="&#92;&#111;&#109;&#101;&#103;&#97;" title="Rendered by QuickLaTeX.com" height="8" width="11" style="vertical-align: 0px;"/>)')
plt.legend()
plt.savefig('magnitude1.eps')
plt.show()
</pre>

Fourier Series

Frequency Response

admin

You might also like

Design Procedure for the Phase Lag Compensator with Example in MATLAB

Policy Iteration Algorithm in Python and Tests with Frozen Lake OpenAI Gym Environment- Reinforcement Learning Tutorial

Interpolation of Function Values in MATLAB