Prześlij pliki do '09/do_sprawdzenia'

This commit is contained in:
Adrianna Hornicka 2021-05-25 23:55:51 +02:00
parent adba0efe20
commit 9f7d02dd2d
2 changed files with 1543 additions and 0 deletions

1345
09/do_sprawdzenia/9_6.ipynb Normal file

File diff suppressed because one or more lines are too long

198
09/do_sprawdzenia/9_7.ipynb Normal file
View File

@ -0,0 +1,198 @@
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Exercise 7.1 \n",
"#### The notebook for this chapter, chap07.ipynb, contains additional examples and explanations. Read through it and run the code."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Done"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Exercise 1\n",
"#### In this chapter, I showed how we can express the DFT and inverse DFT as matrix multiplications. These operations take time proportional to $N^2$, where $N$ is the length of the wave array. That is fast enough for many applications, but there is a faster algorithm, the Fast Fourier Transform (FFT), which takes time proportional to $N \\log N$.\n",
"\n",
"##### The key to the FFT is the Danielson-Lanczos lemma:\n",
"\n",
"$DFT(y)[n] = DFT(e)[n] + \\exp(-2 \\pi i n / N) DFT(o)[n]$\n",
"\n",
"##### Where $ DFT(y)[n]$ is the $n$th element of the DFT of $y$; $e$ is the even elements of $y$, and $o$ is the odd elements of $y$. This lemma suggests a recursive algorithm for the DFT:\n",
"##### Given a wave array, $y$, split it into its even elements, $e$, and its odd elements, $o$.\n",
"##### Compute the DFT of $e$ and $o$ by making recursive calls.\n",
"##### Compute $DFT(y)$ for each value of $n$ using the Danielson-Lanczos lemma.\n",
"##### For the base case of this recursion, you could wait until the length of $y$ is 1. In that case, $DFT(y) = y$. Or if the length of $y$ is sufficiently small, you could compute its DFT by matrix multiplication, possibly using a precomputed matrix.\n",
"##### Hint: I suggest you implement this algorithm incrementally by starting with a version that is not truly recursive. In Step 2, instead of making a recursive call, use dft or np.fft.fft. Get Step 3 working, and confirm that the results are consistent with the other implementations. Then add a base case and confirm that it works. Finally, replace Step 2 with recursive calls.\n",
"##### One more hint: Remember that the DFT is periodic; you might find np.tile useful.\n"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import sys\n",
"sys.setrecursionlimit(10000)"
]
},
{
"cell_type": "code",
"execution_count": 51,
"metadata": {},
"outputs": [],
"source": [
"ys = [-0.2, 0.1, -0.1, -0.5]"
]
},
{
"cell_type": "code",
"execution_count": 53,
"metadata": {},
"outputs": [],
"source": [
"#not truly recursive version\n",
"def FFT(ys):\n",
" N = len(ys)\n",
" if N == 1:\n",
" return ys\n",
" e = np.fft.fft(ys[::2]) #even\n",
" o = np.fft.fft(ys[1::2]) #odd\n",
" n = np.arange(N)\n",
" exp = np.exp(-1j * 2 * np.pi * n / N)\n",
" \n",
" return np.tile(e, 2) + exp * np.tile(o, 2) #Danielson-Lanczos lemma"
]
},
{
"cell_type": "code",
"execution_count": 54,
"metadata": {},
"outputs": [],
"source": [
"#book implementation\n",
"def dft(ys):\n",
" N = len(ys)\n",
" ts = np.arange(N) / N\n",
" fs = np.arange(N)\n",
" args = np.outer(ts, fs)\n",
" M = np.exp(1j * PI2 * args)\n",
" amps = M.conj().transpose().dot(ys)\n",
" return amps"
]
},
{
"cell_type": "code",
"execution_count": 62,
"metadata": {},
"outputs": [],
"source": [
"#resursive\n",
"def FFT_recursive(ys):\n",
" N = len(ys)\n",
" if N == 1:\n",
" return ys \n",
" e = fft(ys[::2])\n",
" o = fft(ys[1::2]) \n",
" n = np.arange(N)\n",
" exp = np.exp(-1j * PI2 * n / N)\n",
" \n",
" return np.tile(e, 2) + exp * np.tile(o, 2)"
]
},
{
"cell_type": "code",
"execution_count": 59,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([-0.7+0.0000000e+00j, -0.1-6.0000000e-01j, 0.1+4.8985872e-17j,\n",
" -0.1+6.0000000e-01j])"
]
},
"execution_count": 59,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"FFT(ys)"
]
},
{
"cell_type": "code",
"execution_count": 60,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([-0.7+0.00000000e+00j, -0.1-6.00000000e-01j, 0.1+1.46957616e-16j,\n",
" -0.1+6.00000000e-01j])"
]
},
"execution_count": 60,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"dft(ys)"
]
},
{
"cell_type": "code",
"execution_count": 63,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([-0.7+0.0000000e+00j, -0.1-6.0000000e-01j, 0.1+4.8985872e-17j,\n",
" -0.1+6.0000000e-01j])"
]
},
"execution_count": 63,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"FFT_recursive(ys)"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.5"
}
},
"nbformat": 4,
"nbformat_minor": 4
}