From 604fa15f94c23eac1797cafa0aacb4d262a57d8b Mon Sep 17 00:00:00 2001 From: pavelkomarov Date: Wed, 15 Jan 2025 16:12:02 -0800 Subject: [PATCH] added dft play, renamed notebooks folder, since it now contains more than just examples, and bumped version number because change of readme to link other places will merit a rerelease --- README.md | 2 +- {examples => notebooks}/chebyshev.ipynb | 0 {examples => notebooks}/dct_types.ipynb | 0 notebooks/dft_play.ipynb | 208 ++++++++++++++++++++++++ {examples => notebooks}/fourier.ipynb | 0 specderiv/__init__.py | 2 +- 6 files changed, 210 insertions(+), 2 deletions(-) rename {examples => notebooks}/chebyshev.ipynb (100%) rename {examples => notebooks}/dct_types.ipynb (100%) create mode 100644 notebooks/dft_play.ipynb rename {examples => notebooks}/fourier.ipynb (100%) diff --git a/README.md b/README.md index 90c338b..929da42 100644 --- a/README.md +++ b/README.md @@ -28,7 +28,7 @@ You should now be able to >>> yth_n = np.sin(th_n) # must be periodic on the domain >>> dyth_n = fourier_deriv(yth_n, 1) ``` -For further usage examples, see the Jupyter notebooks in the examples folder: [Chebyshev](https://github.com/pavelkomarov/spectral-derivatives/blob/main/examples/chebyshev.ipynb) and [Fourier](https://github.com/pavelkomarov/spectral-derivatives/blob/main/examples/fourier.ipynb) +For further usage examples, see the Jupyter notebooks in the notebooks folder: [Chebyshev](https://github.com/pavelkomarov/spectral-derivatives/blob/main/notebooks/chebyshev.ipynb) and [Fourier](https://github.com/pavelkomarov/spectral-derivatives/blob/main/notebooks/fourier.ipynb) ## References diff --git a/examples/chebyshev.ipynb b/notebooks/chebyshev.ipynb similarity index 100% rename from examples/chebyshev.ipynb rename to notebooks/chebyshev.ipynb diff --git a/examples/dct_types.ipynb b/notebooks/dct_types.ipynb similarity index 100% rename from examples/dct_types.ipynb rename to notebooks/dct_types.ipynb diff --git a/notebooks/dft_play.ipynb b/notebooks/dft_play.ipynb new file mode 100644 index 0000000..b9367fe --- /dev/null +++ b/notebooks/dft_play.ipynb @@ -0,0 +1,208 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "738dfead-c1af-4a2e-ba1a-e523cb5190ce", + "metadata": {}, + "source": [ + "# Playing with the DFT\n", + "\n", + "The Discrete Fourier Transform definition can be twisted into several equivalent forms, which may be worth understanding. I use [very similar manipulations to prove $Y_{M-k} = Y_k$ in the math](https://pavelkomarov.com/spectral-derivatives/math.pdf)." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "18531c04-8603-4dc2-b026-63a2b1aef5ea", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "\n", + "def FFT(k):\n", + " \"\"\"just to show DFT(k) definitely == FFT(k) == a_k\"\"\"\n", + " return np.fft.fft(x)[k]\n", + "\n", + "def DFT(k):\n", + "\t\"\"\"standard way to find DFT coefficients\"\"\"\n", + "\ta_k = 0\n", + "\tfor n in range(len(x)):\n", + "\t\ta_k += x[n] * np.exp(-2*np.pi/len(x) * 1j * n * k)\n", + "\treturn a_k\n", + "\n", + "def DFT_(k):\n", + "\t\"\"\"finding DFT coefficients by summing in the opposite direction, using -n\"\"\"\n", + "\ta_k = 0\n", + "\tfor n in range(1, len(x)+1):\n", + "\t\ta_k += x[-n] * np.exp(2*np.pi/len(x) * 1j * n * k)\n", + "\treturn a_k\n", + "\n", + "def DFT__(k, L=3):\n", + "\t\"\"\"finding DFT coefficients with arbitrary roll-around, a mix of the above two formulae\"\"\"\n", + "\ta_k = 0\n", + "\tfor n in range(len(x)-L+1, len(x)+1): # from M-L to M, covers first half\n", + "\t\ta_k += x[-n] * np.exp(2*np.pi/len(x) * 1j * n * k)\n", + "\tfor n in range(L, len(x)):\n", + "\t\ta_k += x[n] * np.exp(-2*np.pi/len(x) * 1j * n * k)\n", + "\treturn a_k\n", + "\n", + "def DFT_even(k):\n", + "\t\"\"\"a formula which equals the other three for *even* signals\"\"\"\n", + "\ta_k = 0\n", + "\tfor n in range(len(x)):\n", + "\t\ta_k += x[n] * np.exp(2*np.pi/len(x) * 1j * n * k)\n", + "\treturn a_k" + ] + }, + { + "cell_type": "markdown", + "id": "d8478cfa-acbc-4573-8e40-4f09e25d8861", + "metadata": {}, + "source": [ + "Let's try it on an even signal. This makes the input vector palindromic: $[x_0, x_1, ... x_N, x_{-(N-1)}, ...x_{-1}]$" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "0aaab8b9-7634-4312-a5e3-0c4cd75bf37f", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "k = 0 , def = 0 : 30.00 - 0.00e+00j\n", + "k = 0 , def = 1 : 30.00 - 0.00e+00j\n", + "k = 0 , def = 2 : 30.00 - 0.00e+00j\n", + "k = 0 , def = 3 : 30.00 - 0.00e+00j\n", + "k = 0 , def = 4 : 30.00 - 0.00e+00j\n", + "k = 1 , def = 0 : -8.00 - 0.00e+00j\n", + "k = 1 , def = 1 : -8.00 + 2.66e-15j\n", + "k = 1 , def = 2 : -8.00 - 1.68e-15j\n", + "k = 1 , def = 3 : -8.00 + 1.78e-15j\n", + "k = 1 , def = 4 : -8.00 - 2.66e-15j\n", + "k = 2 , def = 0 : 6.00 - 0.00e+00j\n", + "k = 2 , def = 1 : 6.00 + 4.44e-15j\n", + "k = 2 , def = 2 : 6.00 - 2.48e-15j\n", + "k = 2 , def = 3 : 6.00 - 4.44e-16j\n", + "k = 2 , def = 4 : 6.00 - 4.44e-15j\n", + "k = 3 , def = 0 : -2.00 - 0.00e+00j\n", + "k = 3 , def = 1 : -2.00 + 7.53e-15j\n", + "k = 3 , def = 2 : -2.00 - 4.59e-15j\n", + "k = 3 , def = 3 : -2.00 + 6.61e-15j\n", + "k = 3 , def = 4 : -2.00 - 7.53e-15j\n", + "k = 4 , def = 0 : 6.00 - 0.00e+00j\n", + "k = 4 , def = 1 : 6.00 + 7.99e-15j\n", + "k = 4 , def = 2 : 6.00 - 4.07e-15j\n", + "k = 4 , def = 3 : 6.00 - 8.88e-16j\n", + "k = 4 , def = 4 : 6.00 - 7.99e-15j\n", + "k = 5 , def = 0 : -8.00 - 0.00e+00j\n", + "k = 5 , def = 1 : -8.00 + 1.47e-14j\n", + "k = 5 , def = 2 : -8.00 - 9.76e-15j\n", + "k = 5 , def = 3 : -8.00 + 1.20e-14j\n", + "k = 5 , def = 4 : -8.00 - 1.47e-14j\n" + ] + } + ], + "source": [ + "x = [4, 3, 5, 10, 5, 3]\n", + "\n", + "def format_complex(z):\n", + "\tpm = \" + \" if z.imag < 0 else \" - \"\n", + "\timag = f\"{np.abs(z.imag):.2e}\" if np.abs(z.imag) < 1e-6 else f\"{np.abs(z.imag):.3f}\"\n", + "\treturn f\"{z.real:.2f}\" + pm + imag + \"j\"\n", + "\n", + "for k in range(0, len(x)):\n", + "\tfor i,f in enumerate([FFT, DFT, DFT_, DFT__, DFT_even]):\n", + "\t\tprint(\"k =\", k, \", def =\", i+1, \":\", format_complex(f(k)))" + ] + }, + { + "cell_type": "markdown", + "id": "4be9a3ba-49db-4f62-824c-f32705ff9e12", + "metadata": {}, + "source": [ + "Note that if we do it with a signal that *isn't* perfectly even, the last definition doesn't equal the others." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "bd707baf-e504-476a-a038-d85bc261ce64", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "k = 0 , def = 1 : 29.00 - 0.00e+00j\n", + "k = 0 , def = 2 : 29.00 - 0.00e+00j\n", + "k = 0 , def = 5 : 29.00 - 0.00e+00j\n", + "k = 1 , def = 1 : -8.50 + 0.866j\n", + "k = 1 , def = 2 : -8.50 + 0.866j\n", + "k = 1 , def = 5 : -8.50 - 0.866j\n", + "k = 2 , def = 1 : 6.50 + 0.866j\n", + "k = 2 , def = 2 : 6.50 + 0.866j\n", + "k = 2 , def = 5 : 6.50 - 0.866j\n", + "k = 3 , def = 1 : -1.00 - 1.11e-16j\n", + "k = 3 , def = 2 : -1.00 + 5.14e-15j\n", + "k = 3 , def = 5 : -1.00 - 5.14e-15j\n", + "k = 4 , def = 1 : 6.50 - 0.866j\n", + "k = 4 , def = 2 : 6.50 - 0.866j\n", + "k = 4 , def = 5 : 6.50 + 0.866j\n", + "k = 5 , def = 1 : -8.50 - 0.866j\n", + "k = 5 , def = 2 : -8.50 - 0.866j\n", + "k = 5 , def = 5 : -8.50 + 0.866j\n" + ] + } + ], + "source": [ + "x = [4, 3, 5, 10, 5, 2]\n", + "\n", + "for k in range(0, len(x)):\n", + "\tfor i,f in zip((1, 2, 5), [FFT, DFT, DFT_even]):\n", + "\t\tprint(\"k =\", k, \", def =\", i, \":\", format_complex(f(k)))" + ] + }, + { + "cell_type": "markdown", + "id": "7cdc97c4-3b5a-4539-a45d-db19d79ff378", + "metadata": {}, + "source": [ + "[\"One, two, five--\" \n", + "\"Three, sir!\" \n", + "\"Three!\"](https://www.youtube.com/watch?v=xOrgLj9lOwk&t=129s)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "84cfeed2-f49c-4b46-a6db-c5d5cf94f8a2", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "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.13.1" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/examples/fourier.ipynb b/notebooks/fourier.ipynb similarity index 100% rename from examples/fourier.ipynb rename to notebooks/fourier.ipynb diff --git a/specderiv/__init__.py b/specderiv/__init__.py index 940de79..886591d 100644 --- a/specderiv/__init__.py +++ b/specderiv/__init__.py @@ -1,3 +1,3 @@ from .specderiv import cheb_deriv, fourier_deriv -__version__ = '0.4' \ No newline at end of file +__version__ = '0.5' \ No newline at end of file