Skip to content

[WIP] add new bonus section regarding "normalization" in Fourier Transform #39

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
NimaSarajpoor wants to merge 3 commits into
base: main
Choose a base branch
from
+212 −0
Open

[WIP] add new bonus section regarding "normalization" in Fourier Transform #39

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
e5a7607
add new bonus section
NimaSarajpoor Jul 26, 2026
433d455
revised bonus section
NimaSarajpoor Jul 26, 2026
618c0a7
minor changes
NimaSarajpoor Jul 30, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
212 changes: 212 additions & 0 deletions docs/Tutorial_FFT_Based_SDP.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@
"outputs": [],
"source": [
"import numpy as np\n",
"import pyfftw\n",
"\n",
"from IPython.display import Image"
]
Expand Down Expand Up @@ -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": {
Expand Down
Loading