"*[CEGM1000 MUDE](http://mude.citg.tudelft.nl/): Week 2.3, Signal Processing. For: November 27, 2024*"
]
},
{
"cell_type": "markdown",
"id": "eff9791a",
"metadata": {},
"source": [
"The goal of this workshop to work with the _Discrete Fourier Transform_ (DFT), implemented in Python as the _Fast Fourier Transform_ (FFT) through `np.fft.fft`, and to understand and interpret its output.\n",
"\n",
"The notebook consists of two parts:\n",
"- The first part (Task 0) is a demonstration of the use of the DFT (_you read and execute the code cells_),\n",
"- The second part is a simple exercise with the DFT (_you write the code_).\n",
"\n",
"To start off, let's do a quick quiz question: _what is the primary purpose of the DFT?_\n",
"That's right! We convert our signal into the frequency domain.\n",
"\n",
"And if you would like an additional explanation of the key frequencies, you can find it [here](https://medium.com/@kovalenko.alx/fun-with-fourier-591662576a77)."
]
},
{
"cell_type": "markdown",
"id": "bfca25a5-75ce-4837-9387-01f95be10bd0",
"metadata": {
"id": "0491cc69"
},
"source": [
"<div style=\"background-color:#facb8e; color: black; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px; width: 95%\"> <p>Note the use of <code>zip</code>, <code>stem</code>, <code>annotate</code> and the modulo operator\n",
"<code>%</code>. Refer to PA 2.3 if you do not understand these tools. Furthermore, note that the term _modulus_ is also used here (and in the textbook), which is another term for _absolute value._</p></div>"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "da13fbf3",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"import numpy as np\n",
"from matplotlib import pyplot as plt"
]
},
{
"cell_type": "markdown",
"id": "13234855",
"metadata": {},
"source": [
"## Task 0: Demonstration of DFT using pulse function\n",
"\n",
"In the first part of this notebook, we use $x(t)=\\Pi(\\frac{t}{4})$, and its Fourier transform $X(f)=4 \\,\\textrm{sinc}(4f)$, as an example (see the first worked example in Chapter 3 on the Fourier transform). The pulse lasts for 4 seconds in the time domain; for convenience, below it is not centered at $t$ = 0, but shifted (delayed) to the right.\n",
"\n",
"The signal $x(t)$ clearly is non-periodic and is an energy signal; apart from a short time span of 'activity' it is zero elsewhere.\n",
"\n",
"We create a pulse function $x(t)$ in discrete time $x_n$ by numpy. The total signal duration is $T$ = 20 seconds (observation or record length). The sampling interval is $\\Delta t$ = 1 second. There are $N$ = 20 samples, and each sample represents 1 second, hence $N \\Delta t = T$.\n",
"\n",
"Note that the time array starts at $t$ = 0 s, and hence, the last sample is at $t$ = 19 s (and not at $t$ = 20 s, as then we would have 21 samples).\n",
"### Task 0.2: Evaluate (and visualize) the DFT\n",
"\n",
"We use `numpy.fft.fft` to compute the one-dimensional Discrete Fourier Transform (DFT) of $x_n$, which takes a signal as argument and returns an array of coefficients (refer to the [documentation](https://numpy.org/doc/stable/reference/generated/numpy.fft.fft.html) as needed).\n",
"\n",
"The DFT converts $N$ samples of the time domain signal $x_n$, into $N$ samples of the frequency domain. In this case, it produces $X_k$ with $k = 0, 1, 2, \\ldots, N-1$, with $N$ = 20, which are complex numbers. "
"Read the code cell below before executing it and identify the following:\n",
"<ol>\n",
" <li>Where is the DFT computed, and what is the output?</li>\n",
" <li>Why is the modulus (absolute value) used on the DFT output?</li>\n",
" <li>What are the values are used for the x and y axes of the plot?</li>\n",
" <li>How is frequency information found and added to the plot (mathematically)?</li>\n",
"</ol>\n",
"Once you understand the figure, continue reading to understand <em>what do these 20 complex numbers mean, and how should we interpret them?</em>\n",
"</p>\n",
"</div>"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "3acb465b",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"abs_fft = np.abs(np.fft.fft(xt))\n",
"index_fft = np.arange(0,20,1)\n",
"plt.plot(index_fft, abs_fft, 'o')\n",
"\n",
"freq = np.arange(0, 1, 0.05)\n",
"for x,y in zip(index_fft, abs_fft):\n",
" if x%5 == 0 or x==19:\n",
" label = f\"f={freq[x]:.2f} Hz\"\n",
" plt.annotate(label, \n",
" (x,y),\n",
" textcoords=\"offset points\", \n",
" xytext=(0,10),\n",
" fontsize=10,\n",
" ha='center') \n",
"plt.ylim(0,5)\n",
"plt.xlim(-2,21)\n",
"plt.xlabel('fft-index')\n",
"plt.ylabel('$|X_k|$')\n",
"plt.stem(index_fft, abs_fft);"
]
},
{
"cell_type": "markdown",
"id": "09c72cdd",
"metadata": {},
"source": [
"The frequency resolution $\\Delta f$ equals one-over-the-measurement-duration, hence $\\Delta f = 1/T$. With that knowledge, we can reconstruct the frequencies expressed in Hertz.\n",
"\n",
"The spectrum of the sampled signal is periodic in the sampling frequency $f_s$ which equals 1 Hz ($\\Delta t$ = 1 s). Therefore, it is computed just for one period $[0,f_s)$. The last value, with index 19, represents the component with frequency $f$ = 0.95 Hz. A spectrum is commonly presented and interpreted as double-sided, so in the above graph we can interpret the spectral components with indices 10 to 19 corresponding to negative frequencies. \n",
"Now we can interpret the DFT of the $x(t)=\\Pi(\\frac{t}{4})$. The sampling interval is $\\Delta t$ = 1 second, and the obsesrvation length is T=20 seconds, and we have N=20 samples. \n",
"- We can recognize a bit of a sinc function.\n",
"- The DFT is computed from 0 Hz, for positive frequencies up to $\\frac{fs}{2} = 0.5$ Hz, after which the negative frequencies follow from -0.5 to -0.05 Hz.\n",
"- $X(f)=4 \\textrm{sinc}(4f)$ has its first null at $f=0.25$Hz.\n",
"\n",
"### Task 0.4: Create symmetric plot\n",
"\n",
"For convenient visualization, we may want to explicitly shift the negative frequencies to the left-hand side to create a symmetric plot. We can use `numpy.fft.fftshift` to do that. In other words, the zero-frequency component appears in the center of the spectrum. Now, it nicely shows a symmetric sprectrum. It well resembles the sinc-function (taking into account that we plot the modulus/absolute value). The output of the DFT still consists of $N$ = 20 elements. To enable this, we set up a new frequency array (in Hz) on the interval $[-fs/2,fs/2)$."
"Read the code cell below before executing it and identify how the plot is modified based on the (new) specification of frequency. Note that it is more than just the <code>freq</code> variable!\n",
"### Task 0.5: Showing spectrum only for positive frequencies\n",
"\n",
"In practice, because of the symmetry, one typically plots only the right-hand side of the (double-sided) spectrum, hence only the part for positive frequencies $f \\geq 0$. This is simply a matter of preference, and a way to save some space."
"Confirm that you understand how we have arrived at the plot above, which illustrates the magnitude (amplitude) spectrum for frequencies $f \\in [0,f_s/2)$, rather than $[0,f_s)$.\n",
"</p>\n",
"</div>"
]
},
{
"cell_type": "markdown",
"id": "9f5558b4",
"metadata": {},
"source": [
"## Task 1: Application of DFT using simple cosine"
]
},
{
"cell_type": "markdown",
"id": "46dd5c0c",
"metadata": {},
"source": [
"It is always a good idea, in spectral analysis, to run a test with a very simple, basic signal. In this way you can test and verify your coding and interpretation of the results. \n",
"\n",
"Our basic signal is just a plain cosine. We take the amplitude equal to one, and zero initial phase, so the signal reads $x(t) = \\cos(2 \\pi f_c t)$, with $f_c$ = 3 Hz in this exercise. With such a simple signal, we know in advance how the spectrum should look like. Namely just a spike at $f$ = 3 Hz, and also one at $f$ = -3 Hz, as we're, for mathematical convenience, working with double sided spectra. The spectrum should be zero at all other frequencies.\n",
"\n",
"As a side note: the cosine is strictly a periodic function, not a-periodic (as above); still, the Fourier transform of the cosine is defined as two Dirac delta pulses or peaks (at 3 Hz and -3 Hz). You may want to check out the second worked example in Chapter 3 on the Fourier transform: Fourier transform in the limit."
"Create a sampled (discrete time) cosine signal by sampling at $f_s$ = 10 Hz, for a duration of $T$ = 2 seconds (make sure you use exactly $N$ = 20 samples). Plot the sampled signal, compute its DFT and plot its magnitude spectrum $|X_k|$ with proper labeling of the axes (just like we did in the first part of this notebook). Include a plot of the spectrum of the sampled cosine signal using only the positive frequencies (up to $f_s/2$, as in the last plot of the previous task).\n",
"\n",
"<em>Note: you are expected to produce three separate plots.</em>\n",
*[CEGM1000 MUDE](http://mude.citg.tudelft.nl/): Week 2.3, Signal Processing. For: November 27, 2024*
%% Cell type:markdown id:eff9791a tags:
The goal of this workshop to work with the _Discrete Fourier Transform_ (DFT), implemented in Python as the _Fast Fourier Transform_ (FFT) through `np.fft.fft`, and to understand and interpret its output.
The notebook consists of two parts:
- The first part (Task 0) is a demonstration of the use of the DFT (_you read and execute the code cells_),
- The second part is a simple exercise with the DFT (_you write the code_).
To start off, let's do a quick quiz question: _what is the primary purpose of the DFT?_
That's right! We convert our signal into the frequency domain.
And if you would like an additional explanation of the key frequencies, you can find it [here](https://medium.com/@kovalenko.alx/fun-with-fourier-591662576a77).
<divstyle="background-color:#facb8e; color: black; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px; width: 95%"><p>Note the use of <code>zip</code>, <code>stem</code>, <code>annotate</code> and the modulo operator
<code>%</code>. Refer to PA 2.3 if you do not understand these tools. Furthermore, note that the term _modulus_ is also used here (and in the textbook), which is another term for _absolute value._</p></div>
%% Cell type:code id:da13fbf3 tags:
``` python
importnumpyasnp
frommatplotlibimportpyplotasplt
```
%% Cell type:markdown id:13234855 tags:
## Task 0: Demonstration of DFT using pulse function
In the first part of this notebook, we use $x(t)=\Pi(\frac{t}{4})$, and its Fourier transform $X(f)=4 \,\textrm{sinc}(4f)$, as an example (see the first worked example in Chapter 3 on the Fourier transform). The pulse lasts for 4 seconds in the time domain; for convenience, below it is not centered at $t$ = 0, but shifted (delayed) to the right.
The signal $x(t)$ clearly is non-periodic and is an energy signal; apart from a short time span of 'activity' it is zero elsewhere.
We create a pulse function $x(t)$ in discrete time $x_n$ by numpy. The total signal duration is $T$ = 20 seconds (observation or record length). The sampling interval is $\Delta t$ = 1 second. There are $N$ = 20 samples, and each sample represents 1 second, hence $N \Delta t = T$.
Note that the time array starts at $t$ = 0 s, and hence, the last sample is at $t$ = 19 s (and not at $t$ = 20 s, as then we would have 21 samples).
We use `numpy.fft.fft` to compute the one-dimensional Discrete Fourier Transform (DFT) of $x_n$, which takes a signal as argument and returns an array of coefficients (refer to the [documentation](https://numpy.org/doc/stable/reference/generated/numpy.fft.fft.html) as needed).
The DFT converts $N$ samples of the time domain signal $x_n$, into $N$ samples of the frequency domain. In this case, it produces $X_k$ with $k = 0, 1, 2, \ldots, N-1$, with $N$ = 20, which are complex numbers.
Read the code cell below before executing it and identify the following:
<ol>
<li>Where is the DFT computed, and what is the output?</li>
<li>Why is the modulus (absolute value) used on the DFT output?</li>
<li>What are the values are used for the x and y axes of the plot?</li>
<li>How is frequency information found and added to the plot (mathematically)?</li>
</ol>
Once you understand the figure, continue reading to understand <em>what do these 20 complex numbers mean, and how should we interpret them?</em>
</p>
</div>
%% Cell type:code id:3acb465b tags:
``` python
abs_fft=np.abs(np.fft.fft(xt))
index_fft=np.arange(0,20,1)
plt.plot(index_fft,abs_fft,'o')
freq=np.arange(0,1,0.05)
forx,yinzip(index_fft,abs_fft):
ifx%5==0orx==19:
label=f"f={freq[x]:.2f} Hz"
plt.annotate(label,
(x,y),
textcoords="offset points",
xytext=(0,10),
fontsize=10,
ha='center')
plt.ylim(0,5)
plt.xlim(-2,21)
plt.xlabel('fft-index')
plt.ylabel('$|X_k|$')
plt.stem(index_fft,abs_fft);
```
%% Cell type:markdown id:09c72cdd tags:
The frequency resolution $\Delta f$ equals one-over-the-measurement-duration, hence $\Delta f = 1/T$. With that knowledge, we can reconstruct the frequencies expressed in Hertz.
The spectrum of the sampled signal is periodic in the sampling frequency $f_s$ which equals 1 Hz ($\Delta t$ = 1 s). Therefore, it is computed just for one period $[0,f_s)$. The last value, with index 19, represents the component with frequency $f$ = 0.95 Hz. A spectrum is commonly presented and interpreted as double-sided, so in the above graph we can interpret the spectral components with indices 10 to 19 corresponding to negative frequencies.
Now we can interpret the DFT of the $x(t)=\Pi(\frac{t}{4})$. The sampling interval is $\Delta t$ = 1 second, and the obsesrvation length is T=20 seconds, and we have N=20 samples.
- We can recognize a bit of a sinc function.
- The DFT is computed from 0 Hz, for positive frequencies up to $\frac{fs}{2} = 0.5$ Hz, after which the negative frequencies follow from -0.5 to -0.05 Hz.
- $X(f)=4 \textrm{sinc}(4f)$ has its first null at $f=0.25$Hz.
### Task 0.4: Create symmetric plot
For convenient visualization, we may want to explicitly shift the negative frequencies to the left-hand side to create a symmetric plot. We can use `numpy.fft.fftshift` to do that. In other words, the zero-frequency component appears in the center of the spectrum. Now, it nicely shows a symmetric sprectrum. It well resembles the sinc-function (taking into account that we plot the modulus/absolute value). The output of the DFT still consists of $N$ = 20 elements. To enable this, we set up a new frequency array (in Hz) on the interval $[-fs/2,fs/2)$.
Read the code cell below before executing it and identify how the plot is modified based on the (new) specification of frequency. Note that it is more than just the <code>freq</code> variable!
### Task 0.5: Showing spectrum only for positive frequencies
In practice, because of the symmetry, one typically plots only the right-hand side of the (double-sided) spectrum, hence only the part for positive frequencies $f \geq 0$. This is simply a matter of preference, and a way to save some space.
Confirm that you understand how we have arrived at the plot above, which illustrates the magnitude (amplitude) spectrum for frequencies $f \in [0,f_s/2)$, rather than $[0,f_s)$.
</p>
</div>
%% Cell type:markdown id:9f5558b4 tags:
## Task 1: Application of DFT using simple cosine
%% Cell type:markdown id:46dd5c0c tags:
It is always a good idea, in spectral analysis, to run a test with a very simple, basic signal. In this way you can test and verify your coding and interpretation of the results.
Our basic signal is just a plain cosine. We take the amplitude equal to one, and zero initial phase, so the signal reads $x(t) = \cos(2 \pi f_c t)$, with $f_c$ = 3 Hz in this exercise. With such a simple signal, we know in advance how the spectrum should look like. Namely just a spike at $f$ = 3 Hz, and also one at $f$ = -3 Hz, as we're, for mathematical convenience, working with double sided spectra. The spectrum should be zero at all other frequencies.
As a side note: the cosine is strictly a periodic function, not a-periodic (as above); still, the Fourier transform of the cosine is defined as two Dirac delta pulses or peaks (at 3 Hz and -3 Hz). You may want to check out the second worked example in Chapter 3 on the Fourier transform: Fourier transform in the limit.
Create a sampled (discrete time) cosine signal by sampling at $f_s$ = 10 Hz, for a duration of $T$ = 2 seconds (make sure you use exactly $N$ = 20 samples). Plot the sampled signal, compute its DFT and plot its magnitude spectrum $|X_k|$ with proper labeling of the axes (just like we did in the first part of this notebook). Include a plot of the spectrum of the sampled cosine signal using only the positive frequencies (up to $f_s/2$, as in the last plot of the previous task).
<em>Note: you are expected to produce three separate plots.</em>
