diff --git a/docs/Tutorial_FFT_Based_SDP.ipynb b/docs/Tutorial_FFT_Based_SDP.ipynb index 3774dc8..3c5f356 100644 --- a/docs/Tutorial_FFT_Based_SDP.ipynb +++ b/docs/Tutorial_FFT_Based_SDP.ipynb @@ -8,6 +8,7 @@ "outputs": [], "source": [ "import numpy as np\n", + "import pyfftw\n", "\n", "from IPython.display import Image" ] @@ -289,6 +290,217 @@ "source": [ "Image(filename='array_merge_annotated.gif', width=800, height=300)" ] + }, + { + "cell_type": "markdown", + "id": "578deec0-6b6f-47f4-82c6-8e5457896269", + "metadata": {}, + "source": [ + "## Bonus: Normalization of forward/ backward Fourier Transforms " + ] + }, + { + "cell_type": "markdown", + "id": "b2c31711-fe53-4ac5-9067-9ec4040eee8c", + "metadata": {}, + "source": [ + "(R)FFT and I(R)FFT are implemented in a way that:\n", + "\n", + "* `IFFT(FFT(x)) = x`\n", + "\n", + "* `IRFFT(RFFT(x)) = x`\n", + "\n", + "By convention, (R)FFT is scaled by `1` (i.e. no scaling), and I(R)FFT is scaled by `1 / len(x)`. Several python libraries like numpy, scipy, and pyfftw follow similar convention, i,e. their forward Fourier Transform has no normalization and their backward Fourier Transform scales the tranformation by `1/len(x)`. This default behaviour results in achieving the correct outcome for circular convolution when the computation is done in the frequency domain. This default behaviour also ensures that `IFFT(FFT(x)) == x` and `IRFFT(RFFT(x)) == x`." + ] + }, + { + "cell_type": "markdown", + "id": "7c5734f9-bba5-49d7-bd08-743bff88541d", + "metadata": {}, + "source": [ + "`pyfftw.FFTW` class provides support for the method `execute`, which performs (forward or backward) Fourier Transform but has no normalization. While one can use this thin wrapper to get a performance boost, they should be careful regarding the normalization part and understand whether that needs to be applied or not. In the following examples, we've used numpy's (i)rfft and pyfftw's (i)rfft to better demonstrate this matter." + ] + }, + { + "cell_type": "markdown", + "id": "096b52ea-fbfc-4c3f-a7e8-63d16848a476", + "metadata": {}, + "source": [ + "### Example 1: Use numpy to check IRFFT(RFFT(x)) == x" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "6ac242c6-cb1c-4a43-8449-b2c994fe9c6a", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "\n", + "\n", + "n_T = 4\n", + "T = np.random.rand(n_T)\n", + "T_ref = T.copy()\n", + "T_comp = np.fft.irfft(np.fft.rfft(T))\n", + "np.testing.assert_allclose(T_comp, T_ref, rtol=1e-5, atol=1e-8)" + ] + }, + { + "cell_type": "markdown", + "id": "b7267b69-0d0b-4d32-b962-2810aa318b2c", + "metadata": {}, + "source": [ + "### Example 2: Use pyfftw to check IRFFT(RFFT(x)) == x" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "93c45f2a-99e7-4ba3-b501-d66f987348aa", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import pyfftw\n", + "\n", + "\n", + "n_T = 4\n", + "real_view = pyfftw.empty_aligned(n_T, dtype=\"float64\")\n", + "complex_view = pyfftw.empty_aligned(n_T//2 + 1, dtype=\"complex128\")\n", + "planning_flag = \"FFTW_ESTIMATE\"\n", + "n_threads = 1\n", + "\n", + "rfft_obj = pyfftw.FFTW(\n", + " input_array=real_view,\n", + " output_array=complex_view,\n", + " direction=\"FFTW_FORWARD\",\n", + " flags=(planning_flag,),\n", + " threads=n_threads,\n", + ")\n", + "\n", + "irfft_obj = pyfftw.FFTW(\n", + " input_array=complex_view,\n", + " output_array=real_view,\n", + " direction=\"FFTW_BACKWARD\",\n", + " flags=(planning_flag,),\n", + " threads=n_threads,\n", + ")\n", + "\n", + "T = np.random.rand(n_T)\n", + "T_ref = T.copy()\n", + "\n", + "\n", + "real_view[:] = T\n", + "rfft_obj()\n", + "R = rfft_obj.output_array\n", + "\n", + "complex_view[:] = R.copy()\n", + "irfft_obj()\n", + "T_comp = irfft_obj.output_array # IRFFT(RFFT(T)) \n", + "\n", + "np.testing.assert_allclose(T_comp, T_ref, rtol=1e-5, atol=1e-8)" + ] + }, + { + "cell_type": "markdown", + "id": "73193969-65b8-40c2-b1c9-9619efc1f882", + "metadata": {}, + "source": [ + "### Example 3: Use pyfftw to check IRFFT(RFFT(x)) == x \n", + "However, this time, we will use the thin wrapper \"execute\" that does Fourier transform only without any scaling." + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "id": "7469ecf6-5187-41f0-a339-f938c373bd3b", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Scaling factor: [0.25 0.25 0.25 0.25]\n" + ] + } + ], + "source": [ + "import numpy as np\n", + "import pyfftw\n", + "\n", + "\n", + "n_T = 4\n", + "real_view = pyfftw.empty_aligned(n_T, dtype=\"float64\")\n", + "complex_view = pyfftw.empty_aligned(n_T//2 + 1, dtype=\"complex128\")\n", + "planning_flag = \"FFTW_ESTIMATE\"\n", + "n_threads = 1\n", + "\n", + "rfft_obj = pyfftw.FFTW(\n", + " input_array=real_view,\n", + " output_array=complex_view,\n", + " direction=\"FFTW_FORWARD\",\n", + " flags=(planning_flag,),\n", + " threads=n_threads,\n", + ")\n", + "\n", + "irfft_obj = pyfftw.FFTW(\n", + " input_array=complex_view,\n", + " output_array=real_view,\n", + " direction=\"FFTW_BACKWARD\",\n", + " flags=(planning_flag,),\n", + " threads=n_threads,\n", + ")\n", + "\n", + "T = np.random.rand(n_T)\n", + "T_ref = T.copy()\n", + "\n", + "real_view[:] = T\n", + "rfft_obj.execute()\n", + "R = rfft_obj.output_array\n", + "\n", + "complex_view[:] = R.copy()\n", + "irfft_obj.execute() # gives un-normalized inverse fourier tranformation\n", + "T_comp = irfft_obj.output_array # IRFFT(RFFT(T)), using unnormalized tranformation\n", + "\n", + "scaling_factor = T_ref / T_comp\n", + "print(\"Scaling factor:\", scaling_factor)" + ] + }, + { + "cell_type": "markdown", + "id": "4f5db96d-6c42-4a97-96ea-f425890b9129", + "metadata": {}, + "source": [ + "**We can now check rfft and irfft seperately:**" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "id": "f90da370-4fb1-42a8-90c2-a3bc5d3b1db0", + "metadata": {}, + "outputs": [], + "source": [ + "np.testing.assert_allclose(R, np.fft.rfft(T), rtol=1e-5, atol=1e-8)\n", + "np.testing.assert_allclose(T_comp * 1/len(T) , np.fft.irfft(R), rtol=1e-5, atol=1e-8)" + ] + }, + { + "cell_type": "markdown", + "id": "4b25f85e-d724-4cfe-a6a3-85c97a8801d8", + "metadata": {}, + "source": [ + "As observed, the thin wrapper `execute` in pyfftw does not noramlize the backward fourier tranform, and the scaling factor is `1/len(T) == 0.25`. Therefore, in circular convolution, users need to scale the input or output arrays by `1/n_T` if they use `pyfftw.execute` for computing the backward fourier transform." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ecb981f3-7453-41ad-8545-6070f1c6911e", + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": {