{
"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",
"# Digital sound and Discrete Fourier analysis\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Example 2.5: Thresholding the DFT coefficients\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",
"from forward_compress_reverse import *\n",
"import matplotlib.pyplot as plt\n",
"from scipy.fftpack import dct, idct"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Since the DFT coefficients represent contributions at given frequencies, they are extremely useful for performing operations on sound, and also for compression. The function `forw_comp_rev_DFT`\n",
"in the module `forward_compress_reverse` \n",
"in the library can manipulate the DFT coefficients of the audio sample file in several ways, and plays the result. \n",
"In particular it can set DFT coefficients below a given threshold to zero. The idea behind doing so is that frequency components with small values\n",
"may contribute little to the perception of the sound, so that they can be discarded. This gives a sound with fewer frequencies, \n",
"more suitable for compression. `forw_comp_rev_DFT` accepts the named parameter `threshold` for this purpose. \n",
"This parameter will also force writing the percentage of the DFT coefficients which are zeroed out to the display.\n",
"\n",
"Let us first listen to the result when `threshold=20`."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"x, fs = forw_comp_rev_DFT(threshold=20)\n",
"play(x, fs)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Then for `threshold=70`."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"x, fs = forw_comp_rev_DFT(threshold=70)\n",
"play(x, fs)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Example 2.6: Quantizing the DFT coefficients\n",
"\n",
"\n",
"\n",
"The thresholding from the previous example affects only frequencies with low contribution. \n",
"This makes it too simple for practical use:\n",
"A more neutral way would be to let each DFT coefficient occupy a certain number of bits. This is also called *quantization*, and is closer to what modern audio standards do. \n",
"This is useful if the number of bits we use is less than what we have at the start.\n",
"`forw_comp_rev_DFT` accepts a named parameter `n`, which secures that a DFT coefficient with bit representation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"...d_2d_1d_0.d_{-1}d_{-2}d_{-3}...\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"is truncated so that the bits $d_{n-1},d_{n-2},...$ are discarded. High values of $n$ thus mean more rounding, \n",
"so that one should hear that the sound degrades with increasing $n$. The following three examples verify this:\n",
"Let us first listen to the result for `n=3`."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"x, fs = forw_comp_rev_DFT(n=3)\n",
"play(x, fs)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Then for `n=5`."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"x, fs = forw_comp_rev_DFT(n=5)\n",
"play(x, fs)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Then for `n=7`."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"x, fs = forw_comp_rev_DFT(n=7)\n",
"play(x, fs)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In practice this quantization procedure is also too simple, since the\n",
"human auditory system is more sensitive to certain frequencies, and should thus allocate a higher number of bits for those. \n",
"Modern audio standards take this into account.\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Exercise 2.21: Plotting the interpolant\n",
"\n",
"\n",
"\n",
"Implement code where you do the following:\n",
"\n",
"* at the top you define the function $f(x)=\\cos^6(x)$, and $M=3$,\n",
"\n",
"* compute the unique interpolant from $V_{M,T}$ (i.e. by taking $N=2M+1$ samples over one period), \n",
" as guaranteed by Proposition 9,\n",
"\n",
"* plot the interpolant against $f$ over one period.\n",
"\n",
"Finally run the code also for $M=4$, $M=5$, and $M=6$. Explain why the plots coincide for $M=6$, but not for $M < 6$. Does increasing $M$ above $M=6$ have any effect on the plots?\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"The code can look as follows."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"from numpy import *\n",
"\n",
"f = lambda t:cos(t)**6\n",
"M = 5\n",
"T = 2*pi\n",
"N = 2*M + 1\n",
"t = linspace(0, T, 100)\n",
"x = f(linspace(0, T - T/float(N), N))\n",
"y = fft.fft(x, axis=0)/N\n",
"s = real(y[0])*ones(len(t))\n",
"for k in range(1,int((N+1)/2)): \n",
" s += 2*real(y[k]*exp(2*pi*1j*k*t/float(T)))\n",
"plt.plot(t, s, 'r', t, f(t), 'g')\n",
"plt.legend(['Interpolant from V_{M,T}','f'])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Exercise 2.25: Execution times for the FFT\n",
"\n",
"\n",
"\n",
"Let us compare execution times for the different methods for computing the DFT.\n",
"\n",
"\n",
"**a)**\n",
"Write code which compares the execution times for an $N$-point DFT for the following three cases: Direct implementation of the DFT \n",
"(as in Example 2.4),\n",
"the FFT implementation used in this chapter, and the built-in `fft`-function. Your code should use the sample audio file `castanets.wav`, \n",
"apply the different DFT implementations to the first $N=2^r$ samples of the file for $r=3$ to $r=15$, store the execution times in a vector, and plot these. You can use \n",
"the function `time()` in the `time` module\n",
"to measure the execution time.\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"The code can look as follows."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"x0, fs = audioread('sounds/castanets.wav')\n",
"\n",
"kvals = arange(3,16)\n",
"slowtime = zeros(len(kvals))\n",
"fasttime = zeros(len(kvals))\n",
"fastesttime = zeros(len(kvals))\n",
"N = 2**kvals\n",
"for k in kvals:\n",
" x = x0[0:2**k]\n",
" x = x.astype(complex)\n",
" \n",
" start = time.time()\n",
" dft_impl(x)\n",
" slowtime[k - kvals[0]] = time.time() - start\n",
" \n",
" start = time.time()\n",
" fft_impl(x, fft_kernel_standard)\n",
" fasttime[k - kvals[0]] = time.time() - start\n",
"\n",
" start = time.time()\n",
" fft.fft(x, axis=0)\n",
" fastesttime[k - kvals[0]] = time.time() - start\n",
"\n",
"# a.\n",
"plt.plot(kvals, slowtime, 'ro-', \\\n",
" kvals, fasttime, 'go-', \\\n",
" kvals, fastesttime, 'bo-')\n",
"plt.grid('on')\n",
"plt.title('time usage of the DFT methods')\n",
"plt.legend(['DFT', 'Standard FFT', 'Built-in FFT'])\n",
"plt.xlabel('log_2 N')\n",
"plt.ylabel('time used [s]')\n",
"plt.show()\n",
" \n",
"plt.figure()\n",
"\n",
"# b.\n",
"plt.loglog(N, slowtime, 'ro-', N, fasttime, 'go-', N, fastesttime, 'bo-')\n",
"plt.axis('equal')\n",
"plt.legend(['DFT', 'Standard FFT', 'Built-in FFT'])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"**b)**\n",
"A problem for large $N$ is that there is such a big difference in the execution times between the two implementations. We can address this by using a loglog-plot. Plot $N$ against execution\n",
"times using the function `loglog`. How should the fact that the number of arithmetic operations are $8N^2$ and $5N\\log_2N$ be reflected in the plot?\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"The two different curves you see should have a derivative approximately equal to one and two, respectively.\n",
"\n",
"\n",
"\n",
"**c)**\n",
"It seems that the built-in FFT is much faster than our own FFT implementation, even though they may use similar algorithms. Try to explain what can be the cause of this.\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"There may be several reasons for this. \n",
"One is that Python code runs slowly when compared to native code, which is much used in the built-in FFT. Also, the built-in `fft` has been subject to much more optimization than we have\n",
"covered here.\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Exercise 2.29: Non-recursive FFT algorithm\n",
"\n",
"\n",
"\n",
"Use the factorization in (2.16) to write a kernel function for a non-recursive FFT implementation.\n",
"In your code, perform the matrix multiplications in equation (2.16) \n",
"from right to left in an (outer) for-loop. For each matrix loop through the different blocks on the diagonal in\n",
"an (inner) for-loop. Make sure you have the right number of blocks on the diagonal, each block being on the form"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{pmatrix} I & D_{N/2^k} \\\\ I & -D_{N/2^k} \\end{pmatrix}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"It may be a good idea to start by implementing multiplication with such a simple matrix first as these are the building blocks in the algorithm. \n",
"Do this so that everything is computed in-place. Also compare the execution times with our original FFT algorithm, as we did in \n",
"Exercise 2.25, and try to explain what you see in this comparison.\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"The library version of the non-recursive FFT is as follows"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" def fft_kernel_nonrec(x, forward):\n",
" \"\"\"\n",
" Compute the FFT of x, using a non-recursive FFT algorithm. \n",
" \n",
" x: a bit-reversed version of the input. Should have only one axis\n",
" forward: Whether the FFT or the IFFT is applied\n",
" \"\"\"\n",
" N = len(x)\n",
" sign = -1\n",
" if not forward:\n",
" sign = 1\n",
" D = exp(sign*2*pi*1j*arange(float(N/2))/N)\n",
" nextN = 1\n",
" while nextN < N:\n",
" k = 0\n",
" while k < N:\n",
" xe, xo = x[k:(k + nextN)], x[(k + nextN):(k + 2*nextN)]\n",
" xo *= D[0::(int(N/(2*nextN)))]\n",
" x[k:(k+2*nextN)] = concatenate([xe + xo, xe - xo])\n",
" k += 2*nextN\n",
" nextN *= 2\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If you add the non-recursive algorithm to the code from Exercise 2.25, you will see that the non-recursive algorithm performs much better. There may be several\n",
"reasons for this. First of all, there are no recursive function calls. Secondly, the values in the matrices $D_{N/2}$ are constructed once and for all with\n",
"the non-recursive algorithm. Code which compares execution times for the original FFT algorithm, our non-recursive implementation, and the split-radix algorithm of the next\n",
"exercise, can look as follows:"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"x0, fs = audioread('sounds/castanets.wav') \n",
"\n",
"kvals = arange(3,16)\n",
"slowtime = zeros(len(kvals))\n",
"fasttime = zeros(len(kvals))\n",
"fastesttime = zeros(len(kvals))\n",
"N = 2**kvals\n",
"for k in kvals:\n",
" x = x0[0:2**k]\n",
" x = x.astype(complex)\n",
" \n",
" start = time.time()\n",
" fft_impl(x, fft_kernel_standard)\n",
" slowtime[k - kvals[0]] = time.time() - start\n",
" \n",
" start = time.time()\n",
" fft_impl(x, fft_kernel_nonrec)\n",
" fasttime[k - kvals[0]] = time.time() - start\n",
"\n",
" start = time.time()\n",
" fft_impl(x, fft_kernel_splitradix)\n",
" fastesttime[k - kvals[0]] = time.time() - start\n",
"\n",
"plt.plot(kvals, slowtime, 'ro-', \\\n",
" kvals, fasttime, 'bo-', \\\n",
" kvals, fastesttime, 'go-')\n",
"plt.grid('on')\n",
"plt.title('time usage of the DFT methods')\n",
"plt.legend(['Standard FFT algorithm', \\\n",
" 'Non-recursive FFT', \\\n",
" 'Split radix FFT'])\n",
"plt.xlabel('log_2 N')\n",
"plt.ylabel('time used [s]')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
""
]
}
],
"metadata": {},
"nbformat": 4,
"nbformat_minor": 2
}