{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"# Exercises from Linear algebra, signal processing, and wavelets. A unified approach.\\\\ Python version\n",
"\n",
" **Øyvind Ryan**\n",
"\n",
"Date: **Oct 24, 2018**\n",
"\n",
"\n",
"\n",
"\n",
"# Sound and Fourier series\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Example 1.1: Listen to different channels\n",
"\n",
""
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import os, sys\n",
"sys.path.append(os.path.join(os.getcwd(), 'python'))"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"%matplotlib inline\n",
"\n",
"from numpy import *\n",
"import time\n",
"from sound import *\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The audio sample file we will use is located in the folder `sounds`:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"x, fs = audioread('sounds/castanets.wav')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"It has two sound channels. In such cases `x` is actually a matrix with two columns, each column representing a sound channel. \n",
"To listen to each channel we can run the following code."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"play(x[:, 0], fs)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"play(x[:, 1], fs)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You may not hear a difference between the two channels. \n",
"There may still be differences, however, which only are notable when the channels are sent to different loudspeakers. \n",
"\n",
"We will in the following apply different operations to sound. We will then mostly apply these operations to the sound channels simultaneously.\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Example 1.2: Playing a sound backwards\n",
"\n",
"\n",
"\n",
"At times a popular game has been to play music backwards to try to find secret messages. \n",
"In the old days of analog music this was not so easy, but with digital sound it is quite simple; we just reverse the samples. \n",
"Thus, if $\\boldsymbol{x}=(x_i)_{i=0}^{N-1}$ are the samples of a digital sound, the samples $\\boldsymbol{y}=(y_i)_{i=0}^{N-1}$ of the reverse sound are"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"y_i=x_{N-i-1}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"When we reverse the sound samples, we have to reverse the elements in both sound channels. \n",
"For our audio sample file this can be performed as follows."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"z = x[::(-1), :]\n",
"play(z, fs)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Example 1.3: Playing pure tones.\n",
"\n",
"\n",
"\n",
"You can also create and listen to sound samples on your own, without reading them from file. \n",
"To create the samples of a pure tone (with only one channel) we can write"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"f = 440\n",
"antsec = 3"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"t = linspace(0, antsec, fs*antsec)\n",
"x = sin(2*pi*f*t)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here `f` is the frequency and `antsec` the length in seconds. \n",
"We can now listen to the samples as follows."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"play(x, fs)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You should hear a pleasant sound with a very distinct frequency. \n",
"Most people can only perceive sounds in the range 20 - 20000 Hz: \n",
"You are encouraged to experiment here to see if your hearing matches this range.\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Example 1.4: The square wave\n",
"\n",
"\n",
"\n",
"There are many other ways in which a function can oscillate regularly. An example is the *square wave*, defined by"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"f_s(t)= \n",
"\\begin{cases} 1, & \\text{if $0\\leq t < T/2$}; \\\\ \n",
" -1, & \\text{if $T/2 \\leq t < T$}.\n",
"\\end{cases}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Given a period $T$, it is $1$ on the first half of each period, and $-1$ on the other.\n",
"The square wave with the same period we used for the pure tone can be plotted as follows:"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"T=1/440.\n",
"t=linspace(0,T,100)\n",
" \n",
"totransl=ones(100)*(-1)**(t>=T/2)\n",
"plt.plot(concatenate([t, t+T, t+2*T, t+3*T, t+4*T]),tile(totransl,5), 'k-')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To listen to the square wave, first create the samples for one period."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"antsec = 3\n",
"samplesperperiod = fs/f # The number of samples for one period\n",
"oneperiod = concatenate([ones((samplesperperiod/2),dtype=float), \\\n",
" -ones((samplesperperiod/2),dtype=float)])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here the function `concatenate` in `numpy` was used to combine two vectors to one larger vector. \n",
"This function will be used more in the remainder of the book. \n",
"Then we repeat one period to obtain a sound with the desired length:"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"x = tile(oneperiod, antsec*f) # Repeat one period\n",
"play(x, fs)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We hear a sound which seems to have the same \"base frequency\" as $\\sin(2\\pi 440 t)$,\n",
"but it is less pleasant to listen to: There seems to be some \"sharp corners\" in the sound. \n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Exercise 1.7: The sum of two pure tones\n",
"\n",
"\n",
"\n",
"\n",
"**a)**\n",
"Consider a sum of two pure tones, $f(t)=a\\sin(2\\pi\\nu_1 t)+b\\sin(2\\pi\\nu_2 t)$. For which values of $a,b,\\nu_1,\\nu_2$ is $f$ periodic? What is the\n",
"period of $f$ when it is periodic?\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"$\\sin(2\\pi\\nu_1 t)$ has period $1/\\nu_1$, while $\\sin(2\\pi\\nu_2 t)$ has period $1/\\nu_2$. \n",
"The period is not unique, however. The first one also has period\n",
"$n/\\nu_1$, and the second also $n/\\nu_2$, for any $n$. The sum is periodic if there exist $n_1,n_2$ \n",
"so that $n_1/\\nu_1=n_2/\\nu_2$, i.e. so that there exists a common period between\n",
"the two. This common period will also be a period of $f$. This amounts to that $\\nu_1/\\nu_2=n_1/n_2$, i.e. that $\\nu_1/\\nu_2$ is\n",
"a rational number.\n",
"\n",
"\n",
"\n",
"**b)**\n",
"Find two constant $a$ and $b$ so that the function $f(t)=a\\sin(2\\pi 440 t)+b\\sin(2\\pi 4400 t)$ resembles the right plot of \n",
"Figure 1.1 as closely as\n",
"possible. Generate the samples of this sound, and listen to it.\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"The important thing to note here is that there are two oscillations present in \n",
"the right part of Figure 1.1: \n",
"One slow oscillation with a higher amplitude, and one\n",
"faster oscillation, with a lower amplitude. We see that there are 10 periods of the smaller oscillation within one period of the larger oscillation, so\n",
"that we should be able to reconstruct the figure by using frequencies where one is 10 times the other, such as 440Hz and 4400Hz. Also, we see from\n",
"the figure that the amplitude of the larger oscillation is close to $1$, and close to $0.3$ for the smaller oscillation. A good choice therefore seems to be\n",
"$a=1,b=0.3$. The code can look this:"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"t = arange(0,3,1/float(fs))\n",
"x = sin(2*pi*440*t) + 0.3*sin(2*pi*4400*t)\n",
"x /= abs(x).max()\n",
"play(x, fs)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Exercise 1.8: Playing with different sample rates\n",
"\n",
"\n",
"\n",
"If we provide another sample rate to the `play` functions, \n",
"the sound card will assume a different time distance between neighboring samples. \n",
"Play and listen to the audio sample file again, but with three different sample rates: $2f_s$, $f_s$, and $f_s/2$, \n",
"where $f_s$ is the sample rate returned by `audioread`.\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"The following code can be used"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"x, fs = audioread('sounds/castanets.wav')"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"play(x, fs)"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"play(x, 2*fs)"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"play(x, fs/2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Exercise 1.9: Playing sound with added noise\n",
"\n",
"\n",
"\n",
"Removing noise from recorded sound can be very challenging, but adding noise is simple. \n",
"There are many kinds of noise, but one kind is easily obtained by adding random numbers to the samples of a sound.\n",
"For this we can use the function\n",
"`random.random` as follows."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" z = x + c*(2*random.random(shape(x))-1)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This adds noise to all channels. \n",
"The function for returning random numbers returns numbers between $0$ and $1$, \n",
"and above we have adjusted these so that they are between $-1$ and $1$ instead, as for other sound which can be played by the computer. \n",
"$c$ is a constant (usually smaller than 1) that dampens the noise. \n",
"\n",
"Write code which adds noise to the audio sample file, and listen to the result for damping constants `c=0.4` and `c=0.1`. \n",
"Remember to scale the sound values after you have added noise, since they may be outside $[-1,1]$.\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"The following code can be used."
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"c = 0.1"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"z = x + c*(2*random.random(shape(x))-1)\n",
"z /= abs(z).max()\n",
"play(z, fs)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Exercise 1.10: Playing the triangle wave\n",
"\n",
"\n",
"\n",
"Repeat what you did in Example 1.4, \n",
"but now for the triangle wave of Example 1.5. \n",
"Start by generating the samples for one period, \n",
"then plot five periods, before you generate the sound over a period of three seconds and play it.\n",
"Verify that you generate the same sound as in Example 1.5.\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"The triangle wave can be plotted as follows:"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"totransl=1-4*abs(t-T/2)/T\n",
"plt.plot(concatenate([t, t+T, t+2*T, t+3*T, t+4*T]),tile(totransl,5), 'k-')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The samples for one period are created as follows."
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"oneperiod = concatenate([linspace(-1, 1, samplesperperiod/2), \\\n",
" linspace(1, -1, samplesperperiod/2)])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Then we repeat one period to obtain a sound with the desired length, and play it as follows."
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"x = tile(oneperiod, antsec*f)\n",
"play(x, fs)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Exercise 1.11: Playing the notes in an octave\n",
"\n",
"\n",
"\n",
"In music theory, an octave is a set of pure tones at frequencies $f_0,...,f_{11},f_{12}$ so that the ratio of neighboring tones are the same, and so that $f_{12}$ is double the frequency of $f_0$, i.e. so that"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\frac{f_1}{f_0}=\\frac{f_2}{f_1}=\\cdots=\\frac{f_{12}}{f_{11}}=2^{1/12}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Make a program which plays all the pure tones in an octave, and listen to it with $f_0=440$Hz.\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"The following code can be used."
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"f_s = 44100\n",
"num_sec = 1\n",
"k = 2**(1/12)\n",
"f = 440\n",
"t = linspace(0, num_sec, f_s*num_sec)\n",
"\n",
"for s in range(13):\n",
" play(sin(2*pi*f*t), f_s)\n",
" time.sleep(num_sec)\n",
" f *= k"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Exercise 1.12: The Karplus-Strong algorithm for making guitar-like sounds\n",
"\n",
"\n",
"\n",
"Given initial values $x_0,...,x_p$, the difference equation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"x_{n+p+1}-\\frac{1}{2}(x_{n+1}+x_n) = 0\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"of order $p+1$ is known to create guitar like sounds. Show that all $x_n$ lie in $[-1,1]$ when the initial values do, and\n",
"write a function"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" karplus_strong(x_init, f_s)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"which takes the initial values $(x_0,x_1,...,x_p)$ as input, and plays the resulting sound for ten seconds with sample rate $f_s$. \n",
"Experiment with randomly generated initial values between $-1$ and $1$, as well as different sample rates.\n",
"What can you say about the frequencies in the resulting sound?\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"The following function can be used"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def karplus_strong(x_init, f_s):\n",
" p=shape(x_init)[0] -1\n",
" num_sec = 10\n",
" num_samples = f_s*num_sec\n",
" z = zeros(num_samples)\n",
" z[0:(p+1)] = x_init\n",
" for k in range(p+1, num_samples):\n",
" z[k] = 0.5*(z[k-p]+z[k-p-1])\n",
" play(z, f_s)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This function can be tested as follows"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"p = 100\n",
"f_s = 44100\n",
"x_init = 2*random.random(p + 1) - 1\n",
"\n",
"karplus_strong(x_init, f_s)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Example 1.13: Fourier coefficients of the square wave\n",
"\n",
"\n",
"\n",
"Let us compute the Fourier coefficients of the square wave, as defined in Example 1.4. \n",
"If we first use equation (1.6) we obtain"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"a_0=\\frac{1}{T}\\int_{0}^{T}f_s(t)dt=\\frac{1}{T}\\int_{0}^{T/2}dt-\\frac{1}{T}\\int_{T/2}^Tdt=0.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Using equation (1.7) we get"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{align*}\n",
"a_n &= \\frac{2}{T}\\int_0^T f_s(t) \\cos(2\\pi nt/T)dt \\\\ \n",
" &= \\frac{2}{T}\\int_0^{T/2} \\cos(2\\pi nt/T)dt-\\frac{2}{T}\\int_{T/2}^T\\cos(2\\pi nt/T)dt \\\\ \n",
" &= \\frac{2}{T}\\left[ \\frac{T}{2\\pi n}\\sin(2\\pi nt/T) \\right]_0^{T/2} - \\frac{2}{T}\\left[ \\frac{T}{2\\pi n}\\sin(2\\pi nt/T) \\right]_{T/2}^T \\\\ \n",
" &= \\frac{2}{T}\\frac{T}{2\\pi n}((\\sin(n\\pi)-\\sin 0) - (\\sin(2n\\pi)-\\sin(n\\pi))=0.\n",
"\\end{align*}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally, using equation (1.8) we obtain"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{align*}\n",
"b_n &= \\frac{2}{T}\\int_0^T f_s(t) \\sin(2\\pi nt/T)dt \\\\ \n",
" &= \\frac{2}{T}\\int_0^{T/2} \\sin(2\\pi nt/T)dt-\\frac{2}{T}\\int_{T/2}^T\\sin(2\\pi nt/T)dt \\\\ \n",
" &= \\frac{2}{T}\\left[ -\\frac{T}{2\\pi n}\\cos(2\\pi nt/T) \\right]_0^{T/2} + \\frac{2}{T}\\left[ \\frac{T}{2\\pi n}\\cos(2\\pi nt/T) \\right]_{T/2}^T \\\\ \n",
" &= \\frac{2}{T}\\frac{T}{2\\pi n}((-\\cos(n\\pi)+\\cos 0) + (\\cos(2n\\pi)-\\cos(n\\pi))) \\\\ \n",
" &= \\frac{2(1-\\cos(n\\pi)}{n\\pi}\\\\ \n",
" & = \\begin{cases} 0, & \\text{if $n$ is even}; \\\\ 4/(n\\pi), & \\text{if $n$ is odd}.\\end{cases}\n",
"\\end{align*}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In other words, only the $b_n$-coefficients with $n$ odd in the Fourier series are nonzero. \n",
"This means that the Fourier series of the square wave is"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\frac{4}{\\pi}\\sin(2\\pi t/T) + \\frac{4}{3\\pi}\\sin(2\\pi 3t/T) + \\frac{4}{5\\pi}\\sin(2\\pi 5t/T) + \\frac{4}{7\\pi}\\sin(2\\pi 7t/T) + \\cdots.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"With $N=20$ there are $10$ terms in the sum. \n",
"The corresponding Fourier series can be plotted over one period with the following code."
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"N = 20\n",
"T = 1/440.\n",
"t = linspace(0, T, 100)\n",
"x = zeros(len(t))\n",
"for k in range(1, N + 1, 2):\n",
" x += (4/(k*pi))*sin(2*pi*k*t/T)\n",
"plt.plot(t, x, 'k-')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To see that the Fourier coefficients converge to 0, let us also plot the first $100$ values of $b_n$:"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"k=arange(1,102,2)\n",
"plt.plot(k,4/(k*pi),'k-')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Even though $f$ oscillates regularly between $-1$ and $1$ with period $T$, the discontinuities mean that it is far from the simple $\\sin(2\\pi t/T)$ which corresponds to a pure tone of frequency $1/T$. \n",
"Clearly $b_1\\sin(2\\pi t/T)$ is the dominant term in the Fourier series. \n",
"This is not surprising since the square wave has the same period as this term, and this explains why we heard the same pitch as the pure tone when we listened to the square wave in Example 1.4. \n",
"Additional terms in the Fourier also contribute, and as we include more of these, we gradually approach the square wave. \n",
"The square wave Fourier coefficients decrease as $1/n$, and this pollution makes the sound less pleasant.\n",
"Note also that, there is a sharp jump in the Fourier series near the discontinuity at $T/2$. It turns out that, even as the number of terms in the Fourier series increases, this sharp jump does not decrease in amplitude. \n",
"This kind of behavior near a discontinuity is also called a *Gibbs phenomenon*. \n",
"\n",
"Let us listen to the Fourier series approximations of the square wave."
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"x = tile(x, antsec/T)\n",
"play(x, fs)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Exercise 1.17: Listening to the Fourier series of the triangle wave\n",
"\n",
"\n",
"\n",
"\n",
"**a)**\n",
"Plot the Fourier series of the triangle wave.\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"The following code can be used."
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"t = linspace(0, T, 100)\n",
"x = zeros(len(t))\n",
"for k in range(1, 20, 2):\n",
" x -= (8/(k*pi)**2)*cos(2*pi*k*t/T)\n",
"plt.plot(t, x, 'k-')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"**b)**\n",
"Write code so that you can listen to the Fourier series of the triangle wave. \n",
"How high must you choose $N$ for the Fourier series to be indistinguishable from the triangle wave itself?\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"The following code can be used."
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"x = tile(x, antsec/T)\n",
"play(x, fs)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Exercise 1.23: Fourier series of polynomials\n",
"\n",
"\n",
"\n",
"Write down difference equations for finding the Fourier coefficients of $f(t)=t^{k}$ from those of $f(t)=t^{k-1}$, and write a program which uses this recursion to compute the Fourier coefficients of $f(t)=t^k$. \n",
"Use the program to verify what you computed in Exercise 1.22.\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"Let us define $a_{n,k},b_{n,k}$ as the Fourier coefficients of $t^k$. When $k>0$ and $n>0$, integration by parts gives us the following difference equations:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{align*}\n",
"a_{n,k} &= \\frac{2}{T} \\int_0^T t^{k}\\cos(2\\pi nt/T)dt \\\\ \n",
" &=\\frac{2}{T}\\left( \\left[ \\frac{T}{2\\pi n} t^k\\sin(2\\pi nt/T) \\right]_0^T - \\frac{kT}{2\\pi n}\\int_0^T t^{k-1}\\sin(2\\pi nt/T)dt\\right) \\\\ \n",
" &=-\\frac{kT}{2\\pi n}b_{n,k-1} \\\\ \n",
"b_{n,k} &= \\frac{2}{T} \\int_0^T t^{k}\\sin(2\\pi nt/T)dt \\\\ \n",
" &=\\frac{2}{T}\\left( \\left[ -\\frac{T}{2\\pi n} t^k\\cos(2\\pi nt/T) \\right]_0^T + \\frac{kT}{2\\pi n}\\int_0^T t^{k-1}\\cos(2\\pi nt/T)dt\\right) \\\\ \n",
" &= -\\frac{T^k}{\\pi n} + \\frac{kT}{2\\pi n}a_{n,k-1}.\n",
"\\end{align*}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"When $n>0$, these can be used to express $a_{n,k},b_{n,k}$ in terms of $a_{n,0},b_{n,0}$, for which we clearly have $a_{n,0}=b_{n,0}=0$. \n",
"For $n=0$ we have that $a_{0,k}=\\frac{T^k}{k+1}$ for all $k$.\n",
"\n",
"The following program computes $a_{n,k},b_{n,k}$ recursively when $n>0$."
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def find_fourier_coeffs(n, k, T):\n",
" ank, bnk = 0, 0\n",
" if n==0:\n",
" ank = T**k/(k+1)\n",
" elif k > 0:\n",
" ankprev, bnkprev = find_fourier_coeffs(n, k-1, T)\n",
" ank = -k*T*bnkprev/(2*pi*n)\n",
" bnk = -T**k/(pi*n) + k*T*ankprev/(2*pi*n)\n",
" return ank, bnk"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Exercise 1.24: Fourier series of a given polynomial\n",
"\n",
"\n",
"\n",
"Use the previous exercise to find the Fourier series for $f(x)=-\\frac{1}{3}x^3+\\frac{1}{2}x^2-\\frac{3}{16}x+1$ on the interval $[0,1]$. \n",
"Plot the $9$th order Fourier series for this function. You should obtain the plots from \n",
"Figure 1.5.\n",
"\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"The following code generates the left plot of Figure 1.5."
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"polycoeffs = matrix([[1, -3/16, 1/2, -1/3]])\n",
"num_terms = 9\n",
"\n",
"ank = matrix(zeros(( shape(polycoeffs)[1], num_terms + 1)))\n",
"bnk = ank.copy()\n",
"for k in range(shape(polycoeffs)[1]):\n",
" for n in range(0, num_terms + 1): \n",
" ank[k, n], bnk[k, n] = find_fourier_coeffs(n, k, 1)\n",
"t = arange(0, 1, 0.01); t = t.reshape((1,len(t))); t = matrix(t)\n",
"n = arange(num_terms + 1); n = n.reshape((len(n), 1)); n = matrix(n)\n",
"vals = polycoeffs*(ank*cos(2*pi*n*t) + bnk*sin(2*pi*n*t))\n",
"\n",
"t_plot = array( t.reshape((1, shape(t)[1])) )\n",
"powervals = matrix(zeros(( shape(polycoeffs)[1], shape(t)[1] )))\n",
"for k in range(shape(polycoeffs)[1]): \n",
" powervals[k,:] = t_plot**k\n",
"\n",
"t_plot = array( t.reshape(( shape(t)[1], 1 )) )\n",
"polyvals = polycoeffs*powervals\n",
"plt.plot(t_plot, polyvals.reshape((shape(polyvals)[1],1)),'r')\n",
"plt.plot(t.reshape((shape(t)[1],1)), vals.reshape((shape(vals)[1],1)), 'b')\n",
"plt.axis([0, 1, 0.975, 1])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The polynomial coefficients are stored in the row vector `polycoeffs` at the top. The code is written so that it will work for any polynomial: \n",
"You just need to replace the vector of polynomial coefficients. The Fourier coefficients \n",
"are stored in the matrices `ank`, `bnk`, using the function implemented in the previous exercise. \n",
"Each row in those matrices correspond to Fourier coefficients of a fixed polynomial.\n",
"We also multiplied with matrices containing cosine and sine values at all instances in time.\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Example 1.25: Complex Fourier coefficients of a part of a pure tone\n",
"\n",
"\n",
"\n",
"Let us consider the pure tone $f(t)=e^{2\\pi it/T_2}$ with period $T_2$, but let us consider it only on the interval $[0,T]$ instead, where $T < T_2$, and extended to a periodic function with period $T$. \n",
"The Fourier coefficients are"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{align*}\n",
"y_n &= \\frac{1}{T}\\int_0^T e^{2\\pi it/T_2}e^{-2\\pi int/T}dt = \\frac{1}{2\\pi iT(1/T_2-n/T)} \\left[ e^{2\\pi it(1/T_2-n/T)}\\right]_0^T \\\\ \n",
"&= \\frac{1}{2\\pi i(T/T_2-n)}\\left( e^{2\\pi iT/T_2} - 1 \\right).\n",
"\\end{align*}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here it is only the term $1/(T/T_2-n)$ which depends on $n$, so that $y_n$ can only be large when $n$ is close $T/T_2$. \n",
"let us plot $|y_n|$ for two different combinations of $T,T_2$. First for $T/T_2=0.5$,"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"n=arange(1,21)\n",
"T=0.5\n",
"T2=1.\n",
"plt.plot( n, abs(( exp(2*pi*1j*T/T2) - 1 )/(2*pi*1j*(T/T2-n))) ,'k-')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"then for $T/T_2=0.9$."
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"T=0.9\n",
"T2=1.\n",
"plt.plot( n, abs(( exp(2*pi*1j*T/T2) - 1 )/(2*pi*1j*(T/T2-n))) ,'k-')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In both plots it is seen that many Fourier coefficients contribute, but this is more visible when $T/T_2=0.5$. When $T/T_2=0.9$, most contribution is seen to be\n",
"in the $y_1$-coefficient. This sounds reasonable, since $f$ then is closest to the pure tone $f(t)=e^{2\\pi it/T}$ of frequency $1/T$ (which in turn has $y_1=1$\n",
"and all other $y_n=0$).\n",
"\n",
"Apart from computing complex Fourier series, there is an important lesson to be learned from this example: In order for a periodic function to be\n",
"approximated by other periodic functions, their period must somehow match.\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Example 1.26: Complex Fourier coefficients when there are different frequencies on different parts\n",
"\n",
"\n",
"\n",
"Most sounds change in content over time. Assume that $f$ is equal to a pure tone of frequency $n_1/T$ on $[0,T/2)$, and equal to a pure\n",
"tone of frequency $n_2/T$ on $[T/2,T)$, i.e."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"f(t) = \\begin{cases} e^{2\\pi in_1t/T} & \\text{on } [0,T_2] \\\\ e^{2\\pi in_2t/T} & \\text{ on} [T_2,T) \\end{cases}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"When $n\\neq n_1,n_2$ we have that"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{align*}\n",
"y_n &= \\frac{1}{T} \\left( \\int_0^{T/2} e^{2\\pi in_1t/T}e^{-2\\pi int/T}dt + \\int_{T/2}^T e^{2\\pi in_2t/T}e^{-2\\pi int/T}dt \\right) \\\\ \n",
"&= \\frac{1}{T}\\left( \\left[ \\frac{T}{2\\pi i(n_1-n)}e^{2\\pi i(n_1-n)t/T} \\right]_0^{T/2} \n",
" + \\left[ \\frac{T}{2\\pi i(n_2-n)}e^{2\\pi i(n_2-n)t/T} \\right]_{T/2}^T \\right) \\\\ \n",
"&= \\frac{e^{\\pi i(n_1-n)}-1}{2\\pi i(n_1-n)} + \\frac{1-e^{\\pi i(n_2-n)}}{2\\pi i(n_2-n)}. \n",
"\\end{align*}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let us restrict to the case when $n_1$ and $n_2$ are both even. We see that"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"y_n = \\begin{cases} \n",
" \\frac{1}{2}+\\frac{1}{\\pi i(n_2-n_1)} & n=n_1,n_2\\\\ \n",
" 0 & n\\text{ even }, n\\neq n_1,n_2\\\\ \n",
" \\frac{n_1-n_2}{\\pi i(n_1-n)(n_2-n)} & n\\text{ odd }\n",
" \\end{cases}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here we have computed the cases $n=n_1$ and $n=n_2$ as above. \n",
"Let us plot $|y_n|$ in two cases: First for $n_1=10$, $n_2=12$,"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"lt=39\n",
" \n",
"n1=10 \n",
"n2=12\n",
"yn=zeros(lt)\n",
"yn[n1-1]=abs(1/2.+1/(pi*1j*(n2-n1)))\n",
"yn[n2-1]=yn[n1-1]\n",
"nvals=arange(1,lt+1,2)\n",
"yn[nvals-1]=abs((n1-n2)/(pi*1j*(n1-nvals)*(n2-nvals)))\n",
"plt.plot(arange(1,lt+1),yn,'k-')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"then for $n_1=2$, $n_2=20$."
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"n1=2\n",
"n2=20\n",
"yn=zeros(lt)\n",
"yn[n1-1]=abs(1/2.+1/(pi*1j*(n2-n1)))\n",
"yn[n2-1]=yn[n1-1]\n",
"nvals=arange(1,lt+1,2)\n",
"yn[nvals-1]=abs((n1-n2)/(pi*1j*(n1-nvals)*(n2-nvals)))\n",
"plt.plot(arange(1,lt+1),yn,'k-')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We see that, when $n_1,n_2$ are close, the Fourier coefficients are close to those of a pure tone of frequency $n/T$ with $n\\approx n_1,n_2$, but that also other\n",
"frequencies contribute. When $n_1,n_2$ are further apart, we see that the Fourier coefficients are like the sum of the two base frequencies. Other frequencies contribute also here.\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"There is an important lesson to be learned from this as well: We should be aware of changes in a sound over time, and it may not be smart to use a\n",
"frequency representation over a large interval when we know that there are simpler frequency representations on the smaller intervals.\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Example 1.36: Periodic and symmetric extension\n",
"\n",
"\n",
"\n",
"Let $f$ be the function with period $T$ defined by $f(t)=2t/T-1$ for $0\\leq t < T$. \n",
"In each period the function increases linearly from $-1$ to $1$. Because $f$ is discontinuous at the boundaries, we would expect the Fourier\n",
"series to converge slowly. The Fourier series is a sine-series since $f$ is antisymmetric, and we can compute $b_n$ as"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{align*}\n",
" b_n &= \\frac{2}{T}\\int_0^T \\frac{2}{T}\\left(t-\\frac{T}{2}\\right)\\sin(2\\pi nt/T)dt = \\frac{4}{T^2}\\int_0^{T} \\left(t-\\frac{T}{2}\\right)\\sin(2\\pi nt/T)dt \\\\ \n",
" &= \\frac{4}{T^2}\\int_0^{T} t\\sin(2\\pi nt/T)dt - \\frac{2}{T}\\int_0^{T} \\sin(2\\pi nt/T)dt = -\\frac{2}{\\pi n},\n",
"\\end{align*}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"so that"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"f_N(t) = - \\sum_{n=1}^N \\frac{2}{n\\pi}\\sin(2\\pi nt/T),\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"which indeed converges slowly to $0$. Let us now instead consider the symmetric extension of $f$. Clearly this is the triangle wave with period $2T$, and the Fourier series of this was"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"(\\breve{f})_N(t) = -\\sum_{n\\leq N\\text{, $n$ odd}} \\frac{8}{n^2\\pi^2}\\cos(2\\pi nt/(2T)).\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The second series clearly converges faster than the first, since its Fourier coefficients are $a_n=-8/(n^2\\pi^2)$ (with $n$ odd), while the Fourier coefficients in the first\n",
"series are $b_n=-2/(n\\pi)$. \n",
"\n",
"If we use $T=1/440$, the symmetric extension has period $1/220$, which gives a triangle wave where the first term in the Fourier series has frequency 220Hz. \n",
"Listening to this we should hear something resembling a 220Hz pure tone, since the first term in the Fourier series is the most dominating in the triangle wave.\n",
"Listening to the periodic extension we should hear a different sound. The first term in the Fourier series has frequency 440Hz, but this drowns a bit in the contribution of the other terms in the\n",
"Fourier series, due to the slow convergence of the Fourier series.\n",
"\n",
"The Fourier series with $7$ terms for the periodic- and symmetric extensions of $f$\n",
"can be plotted as follows"
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"N = 7\n",
"T = 1/440.\n",
"t = linspace(0, 4*T, 400)\n",
"y=zeros(400)\n",
"for n in range(1, N + 1):\n",
" y -= (2/(n*pi))*sin(2*pi*n*t/T)\n",
"plt.plot(t, y, 'k-')"
]
},
{
"cell_type": "code",
"execution_count": 38,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"y=zeros(400)\n",
"for n in range(1, N+1, 2):\n",
" y -= 8*cos(2*pi*n*t/(2*T))/( n**2 * pi**2)\n",
"plt.plot(t, y, 'k-')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"It is clear from the plots that the Fourier series for $f$ is\n",
"not a very good approximation, while for the symmetric extension, we cannot differ between the Fourier series and the function.\n",
"\n",
""
]
}
],
"metadata": {},
"nbformat": 4,
"nbformat_minor": 2
}