From e5a7607f582153a1fd247aea81dda01526de0cbb Mon Sep 17 00:00:00 2001 From: NimaSarajpoor Date: Sat, 25 Jul 2026 21:34:01 -0400 Subject: [PATCH 1/3] add new bonus section --- docs/Tutorial_FFT_Based_SDP.ipynb | 161 +++++++++++++++++++++++++++++- 1 file changed, 160 insertions(+), 1 deletion(-) diff --git a/docs/Tutorial_FFT_Based_SDP.ipynb b/docs/Tutorial_FFT_Based_SDP.ipynb index 3774dc8..d0271e3 100644 --- a/docs/Tutorial_FFT_Based_SDP.ipynb +++ b/docs/Tutorial_FFT_Based_SDP.ipynb @@ -2,12 +2,13 @@ "cells": [ { "cell_type": "code", - "execution_count": 1, + "execution_count": 2, "id": "2ea52f7d-9a78-46f2-80b6-5e782785bca4", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", + "import pyfftw\n", "\n", "from IPython.display import Image" ] @@ -289,6 +290,164 @@ "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": [ + "Previously, it was mentioned that circular convolution can be computed as follows: `IFFT(FFT(T) * FFT(Qr_padded))`. Several python libraries like numpy, scipy, and pyfftw, by default, has no normalization in the forward Fourier Transform (FFT or RFFT). However, their implemention for backward Fourier Transform (IFFT or IRFFT) scales the tranformation by `1/N`. This default behaviour results in achieving the correct outcome for circular convolution when the computation is done in the frequency domain. It also ensures that `IFFT(FFT(x))` or `IRFFT(RFFT(x))` is the same as `x`. " + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "6ac242c6-cb1c-4a43-8449-b2c994fe9c6a", + "metadata": {}, + "outputs": [], + "source": [ + "# Example 1\n", + "# Use numpy to check IRFFT(RFFT(x)) == x\n", + "\n", + "import numpy as np\n", + "\n", + "\n", + "n_T = 4\n", + "T = np.random.rand(n_T)\n", + "R = np.fft.rfft(T)\n", + "IR = np.fft.irfft(R)\n", + "np.testing.assert_allclose(T, IR, rtol=1e-5, atol=1e-8)" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "93c45f2a-99e7-4ba3-b501-d66f987348aa", + "metadata": {}, + "outputs": [], + "source": [ + "# Example 2\n", + "# Use pyfftw to check IRFFT(RFFT(x)) == x\n", + "\n", + "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", + "real_view[:] = T\n", + "rfft_obj()\n", + "R = rfft_obj.output_array\n", + "\n", + "complex_view[:] = R.copy()\n", + "irfft_obj()\n", + "# irfft_obj.output_array contains the IRFFT(RFFT(T))\n", + "\n", + "np.testing.assert_allclose(T, irfft_obj.output_array, rtol=1e-5, atol=1e-8)" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "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": [ + "# Example 3\n", + "# Use pyfftw to check IRFFT(RFFT(x)) == x \n", + "# However, this time, we will use the thin wrapper \"execute\"... \n", + "# ... that does Fourier transform only without any scaling.\n", + "\n", + "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", + "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 fourier tranformation\n", + "\n", + "scaling_factor = T / irfft_obj.output_array\n", + "print(\"Scaling factor:\", scaling_factor)" + ] + }, + { + "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 `0.25`, which is $\\frac{1}{n\\_T}$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ecb981f3-7453-41ad-8545-6070f1c6911e", + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { From 433d455755a4ddf3a89ce631a35898f84b85aedf Mon Sep 17 00:00:00 2001 From: NimaSarajpoor Date: Sat, 25 Jul 2026 22:15:59 -0400 Subject: [PATCH 2/3] revised bonus section --- docs/Tutorial_FFT_Based_SDP.ipynb | 24 +++++++++++++++++++----- 1 file changed, 19 insertions(+), 5 deletions(-) diff --git a/docs/Tutorial_FFT_Based_SDP.ipynb b/docs/Tutorial_FFT_Based_SDP.ipynb index d0271e3..5822ba4 100644 --- a/docs/Tutorial_FFT_Based_SDP.ipynb +++ b/docs/Tutorial_FFT_Based_SDP.ipynb @@ -304,12 +304,26 @@ "id": "b2c31711-fe53-4ac5-9067-9ec4040eee8c", "metadata": {}, "source": [ - "Previously, it was mentioned that circular convolution can be computed as follows: `IFFT(FFT(T) * FFT(Qr_padded))`. Several python libraries like numpy, scipy, and pyfftw, by default, has no normalization in the forward Fourier Transform (FFT or RFFT). However, their implemention for backward Fourier Transform (IFFT or IRFFT) scales the tranformation by `1/N`. This default behaviour results in achieving the correct outcome for circular convolution when the computation is done in the frequency domain. It also ensures that `IFFT(FFT(x))` or `IRFFT(RFFT(x))` is the same as `x`. " + "(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": "code", - "execution_count": 4, + "execution_count": 1, "id": "6ac242c6-cb1c-4a43-8449-b2c994fe9c6a", "metadata": {}, "outputs": [], @@ -329,7 +343,7 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": 2, "id": "93c45f2a-99e7-4ba3-b501-d66f987348aa", "metadata": {}, "outputs": [], @@ -377,7 +391,7 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 3, "id": "7469ecf6-5187-41f0-a339-f938c373bd3b", "metadata": {}, "outputs": [ @@ -438,7 +452,7 @@ "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 `0.25`, which is $\\frac{1}{n\\_T}$." + "As observed, the thin wrapper `execute` in pyfftw does not noramlize the backward fourier tranform, and the scaling factor is `0.25`, which is $\\frac{1}{n\\_T}$. 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." ] }, { From 618c0a76537278b7980ee93b3ee99b8153339503 Mon Sep 17 00:00:00 2001 From: NimaSarajpoor Date: Thu, 30 Jul 2026 00:15:33 -0400 Subject: [PATCH 3/3] minor changes --- docs/Tutorial_FFT_Based_SDP.ipynb | 85 ++++++++++++++++++++++--------- 1 file changed, 62 insertions(+), 23 deletions(-) diff --git a/docs/Tutorial_FFT_Based_SDP.ipynb b/docs/Tutorial_FFT_Based_SDP.ipynb index 5822ba4..3c5f356 100644 --- a/docs/Tutorial_FFT_Based_SDP.ipynb +++ b/docs/Tutorial_FFT_Based_SDP.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": 2, + "execution_count": 1, "id": "2ea52f7d-9a78-46f2-80b6-5e782785bca4", "metadata": {}, "outputs": [], @@ -321,36 +321,46 @@ "`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": 1, + "execution_count": 7, "id": "6ac242c6-cb1c-4a43-8449-b2c994fe9c6a", "metadata": {}, "outputs": [], "source": [ - "# Example 1\n", - "# Use numpy to check IRFFT(RFFT(x)) == x\n", - "\n", "import numpy as np\n", "\n", "\n", "n_T = 4\n", "T = np.random.rand(n_T)\n", - "R = np.fft.rfft(T)\n", - "IR = np.fft.irfft(R)\n", - "np.testing.assert_allclose(T, IR, rtol=1e-5, atol=1e-8)" + "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": 2, + "execution_count": 12, "id": "93c45f2a-99e7-4ba3-b501-d66f987348aa", "metadata": {}, "outputs": [], "source": [ - "# Example 2\n", - "# Use pyfftw to check IRFFT(RFFT(x)) == x\n", - "\n", "import numpy as np\n", "import pyfftw\n", "\n", @@ -378,20 +388,32 @@ ")\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", - "# irfft_obj.output_array contains the IRFFT(RFFT(T))\n", + "T_comp = irfft_obj.output_array # IRFFT(RFFT(T)) \n", "\n", - "np.testing.assert_allclose(T, irfft_obj.output_array, rtol=1e-5, atol=1e-8)" + "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": 3, + "execution_count": 13, "id": "7469ecf6-5187-41f0-a339-f938c373bd3b", "metadata": {}, "outputs": [ @@ -404,11 +426,6 @@ } ], "source": [ - "# Example 3\n", - "# Use pyfftw to check IRFFT(RFFT(x)) == x \n", - "# However, this time, we will use the thin wrapper \"execute\"... \n", - "# ... that does Fourier transform only without any scaling.\n", - "\n", "import numpy as np\n", "import pyfftw\n", "\n", @@ -436,23 +453,45 @@ ")\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 fourier tranformation\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 / irfft_obj.output_array\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 `0.25`, which is $\\frac{1}{n\\_T}$. 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." + "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." ] }, {