I'm trying to figure out how to properly normalize the results of a DFT using FFTW. The FFTW tutorial states that the forward (FFTW_FORWARD) discrete Fourier transform of a 1d complex array X of size n computes an array Y, where
Y_k = \sum\limits_{j=0}^{n-1} X_j e^{-2\pi j k \sqrt{-1}/n}
The backward DFT computes:
Y_k = \sum\limits_{j=0}^{n-1} X_j e^{+2\pi j k \sqrt{-1}/n}
These definitions are the same as for real-to-complex transformations.
Furthermore, the tutorial specifies that "FFTW computes an unnormalized transform, in that there is no coefficient in front of the summation in the DFT. In other words, applying the forward and then the backward transform will multiply the input by n." However, it doesn't specify where exactly this re-scaling needs to be done. I suppose this may be application dependant, but am not sure how to use it properly. This answer states that it should be normalized in the forward direction, but I have my doubts, which I will elaborate.
My goal is to figure out how to properly normalize the FFT results in order to get what I expect. So I did a simple 1D transformation first, where I know what to expect exactly: Using the same convention as FFTW (normalisation factor=1, oscillatory factor=-2*pi for the forward fourier transform), when I transform
1/2 (δ(1 + x) - δ(1 - x))
with δ being the dirac delta function, I expect to get:
integral_(-∞)^∞ (1/2 (δ(1 + x) - δ(1 - x))) e^(-2 π i ω x) dx = i sin(2π ω)
the same holds for when I do an IFFT on i sin(2π ω), only now I need to normalize by dividing by n.
Here is the code I use to demonstrate this behaviour:
program use_fftw
use,intrinsic :: iso_c_binding
implicit none
include 'fftw3.f03'
integer, parameter :: N = 1000
integer, parameter :: dp = kind(1.d0)
real(dp), parameter :: pi = 3.1415926d0
real(dp), parameter :: physical_length = 500
real(dp), parameter :: dx = physical_length/real(N)
real(dp), parameter :: dk = 1.d0 / physical_length
integer :: i, ind1, ind2
! for double precision: use double complex & call dfftw_plan_dft_1d
complex(C_DOUBLE_COMPLEX), allocatable, dimension(:) :: arr_out
real(C_DOUBLE), allocatable, dimension(:) :: arr_in
type(C_PTR) :: plan_forward, plan_backward
allocate(arr_in(1:N))
allocate(arr_out(1:N/2+1))
plan_forward = fftw_plan_dft_r2c_1d(N, arr_in, arr_out, FFTW_ESTIMATE)
plan_backward = fftw_plan_dft_c2r_1d(N, arr_out, arr_in, FFTW_ESTIMATE)
!----------------------
! Setup
!----------------------
! add +1: index = 1 corresponds to x=0
ind1 = int(1.d0/dx)+1 ! index where x=1
ind2 = int((physical_length-1.d0)/dx)+1 ! index where x=-1
arr_in = 0
arr_in(ind1) = -0.5d0
arr_in(ind2) = 0.5d0
!----------------------
! Forward
!----------------------
call fftw_execute_dft_r2c(plan_forward, arr_in, arr_out)
write(*,*) "Verification: Max real part of arr_out:", maxval(real(arr_out))
open(unit=666,file='./fftw_output_norm1d_fft.txt', form='formatted')
do i = 1, N/2+1
write(666, '(2E14.5,x)') (i-1)*dk, aimag(arr_out(i))
enddo
close(666)
write(*,*) "Finished! Written results to fftw_output_norm1d_fft.txt"
!----------------------
! Backward
!----------------------
call fftw_execute_dft_c2r(plan_backward, arr_out, arr_in)
arr_in = arr_in/N
open(unit=666,file='./fftw_output_norm1d_real.txt', form='formatted')
do i = 1, N
write(666, '(2E14.5,x)') (i-1)*dx, arr_in(i)
enddo
close(666)
write(*,*) "Finished! Written results to fftw_output_norm1d_real.txt"
deallocate(arr_in, arr_out)
call fftw_destroy_plan(plan_forward)
call fftw_destroy_plan(plan_backward)
end program use_fftw
And the results, perfectly according to what I'd expect:
So in this case, I only normalized (division by n) when going from Fourier space back to real space and obtained what I wanted.
But I ran into problems when I tried to do the same for multiple dimensions. This time, I'm trying to transform
sqrt(π/2) ((δ(-1 + x) - δ(1 + x)) δ(y) + δ(x) (δ(-1 + y) - δ(1 + y)))
which should give
integral_(-∞)^∞ (sqrt(π/2) ((δ(-1 + x) - δ(1 + x)) δ(y) + δ(x) (δ(-1 + y) - δ(1 + y)))) e^(-2 π i {x, y} {a, b}) d{x, y} = +i sin(a) + i sin(b)
I plot the results for x=0 (k_x = 0, respectively):
which seems completely wrong, both in frequency of the sinus wave and the amplitude.
However, transforming back and normalising by dividing by n^2 gives the expected initial conditions, in both x and y direction. Here is the plot for x=0:
I have no idea what I am doing wrong...
Here is the 2d code:
program use_fftw
use,intrinsic :: iso_c_binding
implicit none
include 'fftw3.f03'
integer, parameter :: N = 1000
integer, parameter :: dp = kind(1.d0)
real(dp), parameter :: pi = 3.1415926d0
real(dp), parameter :: physical_length = 500
real(dp), parameter :: dx = physical_length/real(N)
real(dp), parameter :: dk = 1.d0 / physical_length
integer :: i, ind1, ind2
! for double precision: use double complex & call dfftw_plan_dft_1d
complex(C_DOUBLE_COMPLEX), allocatable, dimension(:,:) :: arr_out
real(C_DOUBLE), allocatable, dimension(:,:) :: arr_in
type(C_PTR) :: plan_forward, plan_backward
allocate(arr_in(1:N, 1:N))
allocate(arr_out(1:N/2+1, 1:N))
plan_forward = fftw_plan_dft_r2c_2d(N, N, arr_in, arr_out, FFTW_ESTIMATE)
plan_backward = fftw_plan_dft_c2r_2d(N, N, arr_out, arr_in, FFTW_ESTIMATE)
!----------------------
! Setup
!----------------------
! add +1: index = 1 corresponds to x=0
ind1 = int(1.d0/dx)+1 ! get index where x = 1
ind2 = int((physical_length-1.d0)/dx)+1 ! get index where x = -1
arr_in = 0
! y=0:
arr_in(ind1, 1) = sqrt(pi/2)
arr_in(ind2, 1) = -sqrt(pi/2)
! x=0:
arr_in(1, ind1) = sqrt(pi/2)
arr_in(1, ind2) = -sqrt(pi/2)
!----------------------
! Forward
!----------------------
call fftw_execute_dft_r2c(plan_forward, arr_in, arr_out)
write(*,*) "Verification: Max real part of arr_out:", maxval(real(arr_out))
open(unit=666,file='./fftw_output_norm2d_fft_x=0.txt', form='formatted')
open(unit=667,file='./fftw_output_norm2d_fft_y=0.txt', form='formatted')
do i = 1, N
write(666, '(2E14.5,x)') (i-1)*dk, aimag(arr_out(1,i))
enddo
do i = 1, N/2+1
write(667, '(2E14.5,x)') (i-1)*dk, aimag(arr_out(i,1))
enddo
close(666)
close(667)
write(*,*) "Finished! Written results to fftw_output_normalisation_fft_x.txt and fftw_output_normalisation_fft_y.txt"
!----------------------
! Backward
!----------------------
call fftw_execute_dft_c2r(plan_backward, arr_out, arr_in)
! Normalisation happens here!
arr_in = arr_in/N**2
open(unit=666,file='./fftw_output_norm2d_real_x=0.txt', form='formatted')
open(unit=667,file='./fftw_output_norm2d_real_y=0.txt', form='formatted')
do i = 1, N
write(666, '(2E14.5,x)') (i-1)*dx, arr_in(1, i)
write(667, '(2E14.5,x)') (i-1)*dx, arr_in(i, 1)
enddo
close(666)
close(667)
write(*,*) "Finished! Written results to fftw_output_norm2d_real_x=0.txt and fftw_output_norm2d_real_y=0.txt"
deallocate(arr_in, arr_out)
call fftw_destroy_plan(plan_forward)
call fftw_destroy_plan( plan_backward)
end program use_fftw
and a python plotting tool:
#!/usr/bin/python3
#====================================
# Plots the results of the FFTW
# example programs.
#====================================
import numpy as np
import matplotlib.pyplot as plt
from sys import argv
from time import sleep
errormessage="""
I require an argument: Which output file to plot.
Usage: ./plot_fftw.py <case>
options for case:
1 fftw_output_norm1d_fft.txt
2 fftw_output_norm1d_real.txt
3 fftw_output_norm2d_fft_x=0.txt
4 fftw_output_norm2d_real_x=0.txt
5 fftw_output_norm2d_fft_y=0.txt
6 fftw_output_norm2d_real_y=0.txt
Please select a case: """
#----------------------
# Hardcoded stuff
#----------------------
file_dict={}
file_dict['1'] = ('fftw_output_norm1d_fft.txt', '1d Fourier transform')
file_dict['2'] = ('fftw_output_norm1d_real.txt', '1d Full circle')
file_dict['3'] = ('fftw_output_norm2d_fft_x=0.txt', '2d Fourier transform, x=0')
file_dict['4'] = ('fftw_output_norm2d_real_x=0.txt', '2d Full circle, x=0')
file_dict['5'] = ('fftw_output_norm2d_fft_y=0.txt', '2d Fourier transform, y=0')
file_dict['6'] = ('fftw_output_norm2d_real_y=0.txt', '2d Full circle, y=0')
#------------------------
# Get case from cmdline
#------------------------
case = ''
def enforce_integer():
global case
while True:
case = input(errormessage)
try:
int(case)
break
except ValueError:
print("\n\n!!! Error: Case must be an integer !!!\n\n")
sleep(2)
if len(argv) != 2:
enforce_integer()
else:
try:
int(argv[1])
case = argv[1]
except ValueError:
enforce_integer()
filename,title=file_dict[case]
#-------------------------------
# Read and plot data
#-------------------------------
k, Pk = np.loadtxt(filename, dtype=float, unpack=True)
fig = plt.figure()
ax = fig.add_subplot(111)
# ax.plot(k, Pk, label='power spectrum')
if case in ['1', '3', '5']:
ax.plot(k, Pk, label='recovered wave', lw=3) # ignore negative k
x = np.linspace(k.min(), k.max(), 1000)
if case=='1':
ax.plot(x, np.sin(2*np.pi*x), ':', label='expected wave', lw=3)
if case in ['3', '5']:
ax.plot(x, np.sin(x), ':', label='expected wave', lw=3)
ax.set_title(title)
ax.set_xlabel("k")
ax.set_ylabel("F(k)")
if case in ['2', '4', '6']:
# in this case: k=x, Pk=f(x)
ax.plot(k, Pk, label='recovered original', lw=3) # ignore negative k
N=1000
plen=500
dx=plen/N
x = np.linspace(k.min(), k.max(), 1000)
y = np.zeros(1000)
ind = int(1.0/dx)
if case=='2':
y[ind] = -0.5
y[-ind] = 0.5
if case in ['4', '6']:
y[ind] = np.sqrt(np.pi/2)
y[-ind] = -np.sqrt(np.pi/2)
ax.plot(x, y, ':', label='expected original', lw=3)
ax.set_title(title)
ax.set_xlabel("x")
ax.set_ylabel("f(x)")
ax.legend()
plt.show()
