|
| 1 | +{ |
| 2 | + "cells": [ |
| 3 | + { |
| 4 | + "cell_type": "markdown", |
| 5 | + "metadata": {}, |
| 6 | + "source": [ |
| 7 | + "# The Discrete Fourier Transform\n", |
| 8 | + "\n", |
| 9 | + "## Starting point: the continuous-time case\n", |
| 10 | + "\n", |
| 11 | + "As you recall, a continuous-time periodic function $f(t)$ with angular frequency $\\omega_0=2 \\pi f = 2 \\pi / T$ can be written as \n", |
| 12 | + "\n", |
| 13 | + "$$f(t) = \\sum_{n=-\\infty}^\\infty a_n e^{jn\\omega_0 t}$$\n", |
| 14 | + "\n", |
| 15 | + "Such an expression can be considered as an inverse Fourier Transform (FT), i.e. an operation going from a representation by coefficients $a_n$ in the frequency domain to a representation $f(t)$ in the time-domain.\n", |
| 16 | + "\n", |
| 17 | + "The coefficients $a_n$ are given by\n", |
| 18 | + "\n", |
| 19 | + "$$ a_n = \\frac{1}{T} \\int_0^T f(t) e ^{-jn\\omega_0 t} dt$$ \n", |
| 20 | + "\n", |
| 21 | + "The expression above can be considered as a forward Fourier Transform, i.e. from the time domain to the frequency domain.\n", |
| 22 | + "\n", |
| 23 | + "## Moving to discrete time\n", |
| 24 | + "\n", |
| 25 | + "Now suppose that instead of a continuous-time function $f(t)$, we have a discrete-time signal given e.g. by a set of samples $x[k]=f(k \\Delta t)$, with $\\Delta t$ the sampling period. In that case, we use the Discrete Fourier Transform (DFT) rather than the regular FT. The DFT is defined as:\n", |
| 26 | + "\n", |
| 27 | + "$$X[n] = \\sum_{k=0}^{N-1} x[k] e^{-j 2 \\pi k n /N}$$\n", |
| 28 | + "\n", |
| 29 | + "To better illustrate the relationship between the DFT and FT, let's look at $-j n \\omega_0 t$ in the formula for the FT. This becomes $-j n \\frac{2 \\pi}{T} t = -j n \\frac{2 \\pi}{T} k \\Delta t$. If we need $N$ samples to cover our period $T$, i.e. if $T=N \\Delta t$, this becomes $-j n \\frac{2 \\pi}{N \\Delta t} k \\Delta t = -j n \\frac{2 \\pi}{N} k $. This is exactly what we have in the exponential in the DFT. (Note that this is a rather handwaving argument. The full derivation is more involved and relies on convolutions, delta functions and sync functions.)\n", |
| 30 | + "\n", |
| 31 | + "Also, an expression like $-j 2 \\pi f t = -j 2 \\pi \\frac{n}{N \\Delta t} k \\Delta t$ makes it crystal clear that $k$ is a temporal index, with samples separated in time by $\\Delta t$, whereas $n$ is a frequency index, with different frequencies separated by $ \\frac{1}{N \\Delta t}$.\n", |
| 32 | + "\n", |
| 33 | + "For the inverse DFT, we have\n", |
| 34 | + "\n", |
| 35 | + "$$x[k]=\\frac{1}{N} \\sum_{n=0}^{N-1} X[n] e^{j 2 \\pi k n /N}$$\n", |
| 36 | + "\n", |
| 37 | + "So, our discrete-time signal $x[k]$ is the sum over many frequency contributions. Each term in the sum is a complex exponential $e^{j 2 \\pi k n /N}$ with a complex-valued amplitude $X[n]$. This is very remininscent of what we had with phasors in the continuous time domain. There, we associated a complex amplitude to an oscillation described by a complex exponential $e^{j \\omega t}$. Here, we do the same thing, but in discrete time.\n", |
| 38 | + "\n", |
| 39 | + "(Note: There also exists something called the Discrete-Time Fourier Transform (DTFT), which is similar but different. It is used for aperiodic functions, with intergration limits from $-\\infty$ to $\\infty$, rather than from 0 to $T$.)\n", |
| 40 | + "\n", |
| 41 | + "In case $N$ is a power of 2, there is a very efficient implementation of the DTF, which is called the Fast Fourier Transform (FFT).\n", |
| 42 | + "\n", |
| 43 | + "## First numerical experiments\n", |
| 44 | + "\n", |
| 45 | + "Let's now play around with some numerical experiments of the DFT in NumPy. Be aware that other definitions/conventions in terms of the sign in the exponential or normalisation factors are possible. The conventions used here are typical for engineering. Look up the definition used by NumPy for its DFT implementation [here](https://numpy.org/doc/stable/reference/routines.fft.html#implementation-details). What convention do they use? Note that the name FFT is so ubiquitous that the function called ```fft``` can also be used if $N$ is not a power of 2. It will then fall back on slower algorithms.\n", |
| 46 | + "\n", |
| 47 | + "Now, start from a sine wave with frequency 2 Hz, and construct a discrete version of it by sampling at 10 points per period. Plot one period of the original function, as well as the sampled points." |
| 48 | + ] |
| 49 | + }, |
| 50 | + { |
| 51 | + "cell_type": "code", |
| 52 | + "execution_count": 1, |
| 53 | + "metadata": {}, |
| 54 | + "outputs": [], |
| 55 | + "source": [ |
| 56 | + "# Your code here." |
| 57 | + ] |
| 58 | + }, |
| 59 | + { |
| 60 | + "cell_type": "markdown", |
| 61 | + "metadata": {}, |
| 62 | + "source": [ |
| 63 | + "Based on the formulas above, calculate the frequencies that correspond to the different coefficients in the DFT. Then, calculate the DFT and plot the absolute values of the coefficients. Does this correspond with what you expect?" |
| 64 | + ] |
| 65 | + }, |
| 66 | + { |
| 67 | + "cell_type": "code", |
| 68 | + "execution_count": 3, |
| 69 | + "metadata": {}, |
| 70 | + "outputs": [], |
| 71 | + "source": [ |
| 72 | + "# Your code here." |
| 73 | + ] |
| 74 | + }, |
| 75 | + { |
| 76 | + "cell_type": "markdown", |
| 77 | + "metadata": {}, |
| 78 | + "source": [ |
| 79 | + "This is perhaps a bit bizarre, as we would expect a contribution from both a positive and a negative exponential at frequency 2 Hz, since $\\sin(x)=\\frac{\\exp(j \\omega t) - \\exp(-j \\omega t)}{2j}$. Take a moment to try and figure out what is going on here.\n", |
| 80 | + "\n", |
| 81 | + "Our complex exponential for the $n$-th frequency is $e^{j 2 \\pi k n /N}$, with $k$ and $n$ being integer indices that represent time and frequency respectively. We know that if we add an integer multiple of $2 \\pi$ to the argument of a complex exponential, nothing changes. It's easy to check that replacing $k$ by $k + N$, or $n$ by $n + N$ does exactly that. This results in the time being shifted from $t$ to $t + N \\Delta t = t + T$, and the frequency being shifted from $f$ to $f + N \\frac{1}{N \\Delta t} = f + \\frac{1}{ \\Delta t} $. So, the function is periodic in time with period $T$ (we already knew that!), but it's also periodic in frequency with period $\\frac{1}{ \\Delta t}$. Because this frequency period is the inverse of the sampling time, it's also called the **sampling frequency $f_s$**.\n", |
| 82 | + "\n", |
| 83 | + "Write code to replace the upper half of the frequencies in the example above by equivalent frequencies that are shifted down by $\\frac{1}{\\Delta t}$." |
| 84 | + ] |
| 85 | + }, |
| 86 | + { |
| 87 | + "cell_type": "code", |
| 88 | + "execution_count": 5, |
| 89 | + "metadata": {}, |
| 90 | + "outputs": [], |
| 91 | + "source": [ |
| 92 | + "# You code here." |
| 93 | + ] |
| 94 | + }, |
| 95 | + { |
| 96 | + "cell_type": "markdown", |
| 97 | + "metadata": {}, |
| 98 | + "source": [ |
| 99 | + "There is actually a function ```np.fft.fftfreq``` that generates this array directly. Look up its [documentation](https://numpy.org/doc/stable/reference/generated/numpy.fft.fftfreq.html) and verify it gives the same result as your own code above." |
| 100 | + ] |
| 101 | + }, |
| 102 | + { |
| 103 | + "cell_type": "code", |
| 104 | + "execution_count": 7, |
| 105 | + "metadata": {}, |
| 106 | + "outputs": [], |
| 107 | + "source": [ |
| 108 | + "# Your code here." |
| 109 | + ] |
| 110 | + }, |
| 111 | + { |
| 112 | + "cell_type": "markdown", |
| 113 | + "metadata": {}, |
| 114 | + "source": [ |
| 115 | + "Finally, it would be nice to reorder both the array of the frequencies and the array of FFT coefficients such that zero is in the middle. This can be accomplished by a cyclic shift (i.e. the elements that 'fall out' of the array on the right-hand side, reenter the array from the left), and this is implemented by ```np.fft.fftshift``` ([documentation](https://numpy.org/doc/stable/reference/generated/numpy.fft.fftshift.html)). Do that, and replot the results." |
| 116 | + ] |
| 117 | + }, |
| 118 | + { |
| 119 | + "cell_type": "code", |
| 120 | + "execution_count": 9, |
| 121 | + "metadata": {}, |
| 122 | + "outputs": [], |
| 123 | + "source": [ |
| 124 | + "# Your code here." |
| 125 | + ] |
| 126 | + }, |
| 127 | + { |
| 128 | + "cell_type": "markdown", |
| 129 | + "metadata": {}, |
| 130 | + "source": [ |
| 131 | + "This results in what we intuitively expected." |
| 132 | + ] |
| 133 | + }, |
| 134 | + { |
| 135 | + "cell_type": "markdown", |
| 136 | + "metadata": {}, |
| 137 | + "source": [ |
| 138 | + "# Digging a bit deeper\n", |
| 139 | + "\n", |
| 140 | + "So far, we've only looked at the magnitude of the FFT, but let's now plot the real part and the imaginary part separately. What do you expect?" |
| 141 | + ] |
| 142 | + }, |
| 143 | + { |
| 144 | + "cell_type": "code", |
| 145 | + "execution_count": 11, |
| 146 | + "metadata": {}, |
| 147 | + "outputs": [], |
| 148 | + "source": [ |
| 149 | + "# Your code here." |
| 150 | + ] |
| 151 | + }, |
| 152 | + { |
| 153 | + "cell_type": "markdown", |
| 154 | + "metadata": {}, |
| 155 | + "source": [ |
| 156 | + "Since $\\sin(x)=\\frac{\\exp(j \\omega t) - \\exp(-j \\omega t)}{2j}$, we indeed expect that the coefficients of the complex exponentials should be purely imaginary, which a negative imaginary part of the coefficients for positive frequencies and a positive one for negative frequencies. As for the size of the amplitudes, we expected 0.5 and not 5, but remember that in the formula to go back to the time domain, there is a factor $1/N$, which is 1/10 in this case:\n", |
| 157 | + "\n", |
| 158 | + "$$x[k]=\\frac{1}{N} \\sum_{n=0}^{N-1} X[n] e^{j 2 \\pi k n /N}$$" |
| 159 | + ] |
| 160 | + }, |
| 161 | + { |
| 162 | + "cell_type": "markdown", |
| 163 | + "metadata": {}, |
| 164 | + "source": [ |
| 165 | + "Finally, in case your original time sequence only consists of real values, what can you then conclude in general about the real and imaginary part of the FFT coefficients?\n", |
| 166 | + "\n", |
| 167 | + "Looking back at the formula for the DFT, we have\n", |
| 168 | + "\n", |
| 169 | + "$$X[n] = \\sum_{k=0}^{N-1} x[k] e^{-j 2 \\pi k n /N}$$\n", |
| 170 | + "\n", |
| 171 | + "It's easy to see that $X[-n] = X[n]^*$ in case $x[k]$ is real.\n", |
| 172 | + "\n", |
| 173 | + "This means that $\\Re X[-n] + j \\Im X[-n] = \\Re X[n] - j \\Im X[n]$. So, the real part is even, and the imaginary part is odd, which is what we had in our concrete example above as well." |
| 174 | + ] |
| 175 | + }, |
| 176 | + { |
| 177 | + "cell_type": "markdown", |
| 178 | + "metadata": {}, |
| 179 | + "source": [ |
| 180 | + "## Playing around with frequency range and resolution\n", |
| 181 | + "\n", |
| 182 | + "First, rewrite the code to plot the samples and the magnitude of the FFT as a function, so that you can more easily use it to explore different scenarios. You can use e.g. the following signature:" |
| 183 | + ] |
| 184 | + }, |
| 185 | + { |
| 186 | + "cell_type": "code", |
| 187 | + "execution_count": 13, |
| 188 | + "metadata": {}, |
| 189 | + "outputs": [], |
| 190 | + "source": [ |
| 191 | + "def plot_dft(f_t, T, N):\n", |
| 192 | + " \n", |
| 193 | + " # Your code to plot the sampled function f_t.\n", |
| 194 | + " \n", |
| 195 | + " # Your code to to plot the DFT.\n", |
| 196 | + " \n", |
| 197 | + " pass\n", |
| 198 | + "\n", |
| 199 | + "f = 2\n", |
| 200 | + "plot_dft(lambda t : np.sin(2*np.pi*f*t), T=0.5, N=10)" |
| 201 | + ] |
| 202 | + }, |
| 203 | + { |
| 204 | + "cell_type": "markdown", |
| 205 | + "metadata": {}, |
| 206 | + "source": [ |
| 207 | + "Starting from the example we have studied above, investigate what happens if you sample two periods of the function rather than just one, but keeping the same distance between samples. What do you expect based on the formulas? Do the numerical experiments confirm this?" |
| 208 | + ] |
| 209 | + }, |
| 210 | + { |
| 211 | + "cell_type": "code", |
| 212 | + "execution_count": 15, |
| 213 | + "metadata": {}, |
| 214 | + "outputs": [], |
| 215 | + "source": [ |
| 216 | + "# Your code here." |
| 217 | + ] |
| 218 | + }, |
| 219 | + { |
| 220 | + "cell_type": "markdown", |
| 221 | + "metadata": {}, |
| 222 | + "source": [ |
| 223 | + "From our discussion above, we know that the frequencies are $\\frac{n}{N \\Delta t}$, and have a range $\\frac{1}{\\Delta t}$ after which they repeat. \n", |
| 224 | + "\n", |
| 225 | + "In this case, $\\Delta t$ stays the same, and the plots indeed confirm that the range of frequencies hasn't changed. However, $N$ has increased, which results in frequencies that are closer together. I.e., the frequency resolution has increased." |
| 226 | + ] |
| 227 | + }, |
| 228 | + { |
| 229 | + "cell_type": "markdown", |
| 230 | + "metadata": {}, |
| 231 | + "source": [ |
| 232 | + "Next up, restrict your sampling back to a single period, but this time halve the sampling interval. What happens?" |
| 233 | + ] |
| 234 | + }, |
| 235 | + { |
| 236 | + "cell_type": "code", |
| 237 | + "execution_count": 17, |
| 238 | + "metadata": {}, |
| 239 | + "outputs": [], |
| 240 | + "source": [ |
| 241 | + "# Your code here." |
| 242 | + ] |
| 243 | + }, |
| 244 | + { |
| 245 | + "cell_type": "markdown", |
| 246 | + "metadata": {}, |
| 247 | + "source": [ |
| 248 | + "Since $\\Delta t$ is halved, the range of frequencies is doubled. At the same time $N$ is doubled, so these two effects cancel in the formula for the resolution, which therefore stays the same.\n", |
| 249 | + "\n", |
| 250 | + "So, **to summarise the important conclusions** we've seen:\n", |
| 251 | + "\n", |
| 252 | + "- The frequencies in the DFT periodically repeat with a period equal to the sampling frequency $f_s=1/\\Delta t$. If we choose to have zero as the center frequency, this results in a frequency range of $[-f_s/2, f_s/2]$.\n", |
| 253 | + "- The frequency resolution is given by $\\Delta f = \\frac{1}{N \\Delta t} = \\frac{1}{N \\frac{T}{N}} = \\frac{1}{T}$." |
| 254 | + ] |
| 255 | + }, |
| 256 | + { |
| 257 | + "cell_type": "markdown", |
| 258 | + "metadata": {}, |
| 259 | + "source": [ |
| 260 | + "## Aliasing and the Nyquist criterion" |
| 261 | + ] |
| 262 | + }, |
| 263 | + { |
| 264 | + "cell_type": "markdown", |
| 265 | + "metadata": {}, |
| 266 | + "source": [ |
| 267 | + "We started out by doing the DFT of a sine wave with frequency $f=2$ Hz, using $N=10$ samples to cover a duration of $T=1/f=0.5$ seconds. Since the time interval $\\Delta t$ between samples is $T/N=0.05$ seconds, we know that the frequencies in the DFT are only defined up to multiples of $1/\\Delta t = 20$ Hz.\n", |
| 268 | + "\n", |
| 269 | + "Let's now use the same numerical parameters $N$ and $T$, and do the DFT of a sine wave of frequency 22 Hz. (Note that we'll still use $T=0.5$ seconds, so $T$ is not related to the frequency of the new sine wave anymore.)\n", |
| 270 | + "\n", |
| 271 | + "Update the parameters of your function to plot the results, and try to explain them." |
| 272 | + ] |
| 273 | + }, |
| 274 | + { |
| 275 | + "cell_type": "code", |
| 276 | + "execution_count": 19, |
| 277 | + "metadata": {}, |
| 278 | + "outputs": [], |
| 279 | + "source": [ |
| 280 | + "# Your code here" |
| 281 | + ] |
| 282 | + }, |
| 283 | + { |
| 284 | + "cell_type": "markdown", |
| 285 | + "metadata": {}, |
| 286 | + "source": [ |
| 287 | + "So, rather than ending up with a frequency of 22 Hz, which would be beyond our frequency range of 20 Hz anyhow, we end up with a frequency of 22-20=2 Hz. The time domain graph also clearly shows that the samples lie on a sine wave with a much lower frequency than the frequency of the original signal.\n", |
| 288 | + "\n", |
| 289 | + "Let's now try a freqency of 12 Hz. 12 is smaller than 20, right? So, hopefully things go better now..." |
| 290 | + ] |
| 291 | + }, |
| 292 | + { |
| 293 | + "cell_type": "code", |
| 294 | + "execution_count": null, |
| 295 | + "metadata": {}, |
| 296 | + "outputs": [], |
| 297 | + "source": [ |
| 298 | + "# Your code here" |
| 299 | + ] |
| 300 | + }, |
| 301 | + { |
| 302 | + "cell_type": "markdown", |
| 303 | + "metadata": {}, |
| 304 | + "source": [ |
| 305 | + "Rather than 12 Hz, we end up with 8 Hz. However, if we do 12-20, we get -8 Hz, so this mix-up is still consistent with having only a limited frequency range of 20 Hz that periodically repeats.\n", |
| 306 | + "\n", |
| 307 | + "The fact that high frequencies end up being mistaken for lower frequencies is called **aliasing**.\n", |
| 308 | + "\n", |
| 309 | + "Based on the examples we've seen so far, it seems that because of a limited range of frequencies $[-f_s/2,f_s/2]$ that periodically repeats, in order to get reliable results, the frequency of the signal that we're interested in should be lower than $f_s/2$.\n", |
| 310 | + "\n", |
| 311 | + "What if the signal frequency is exactly $f_s/2$? " |
| 312 | + ] |
| 313 | + }, |
| 314 | + { |
| 315 | + "cell_type": "code", |
| 316 | + "execution_count": null, |
| 317 | + "metadata": {}, |
| 318 | + "outputs": [], |
| 319 | + "source": [ |
| 320 | + "# Your code here." |
| 321 | + ] |
| 322 | + }, |
| 323 | + { |
| 324 | + "cell_type": "markdown", |
| 325 | + "metadata": {}, |
| 326 | + "source": [ |
| 327 | + "This is a pathological edge case: the samples all end up at zero.\n", |
| 328 | + "\n", |
| 329 | + "So, to summarise: **the DFT only gives reliable, i.e. non-aliased results for signal frequencies $f \\lt f_s/2$**.\n", |
| 330 | + "\n", |
| 331 | + "This is called the **Nyquist criterion**, and $f_s/2$ is the Nyquist frequency. Rather than Nyquist, you sometimes can also hear the name of **Shannon** instead. " |
| 332 | + ] |
| 333 | + } |
| 334 | + ], |
| 335 | + "metadata": { |
| 336 | + "kernelspec": { |
| 337 | + "display_name": "Python 3 (ipykernel)", |
| 338 | + "language": "python", |
| 339 | + "name": "python3" |
| 340 | + }, |
| 341 | + "language_info": { |
| 342 | + "codemirror_mode": { |
| 343 | + "name": "ipython", |
| 344 | + "version": 3 |
| 345 | + }, |
| 346 | + "file_extension": ".py", |
| 347 | + "mimetype": "text/x-python", |
| 348 | + "name": "python", |
| 349 | + "nbconvert_exporter": "python", |
| 350 | + "pygments_lexer": "ipython3", |
| 351 | + "version": "3.11.0" |
| 352 | + }, |
| 353 | + "vscode": { |
| 354 | + "interpreter": { |
| 355 | + "hash": "5238573367df39f7286bb46f9ff5f08f63a01a80960060ce41e3c79b190280fa" |
| 356 | + } |
| 357 | + } |
| 358 | + }, |
| 359 | + "nbformat": 4, |
| 360 | + "nbformat_minor": 2 |
| 361 | +} |
0 commit comments