How to normalize the spectrum of a numpy (real) fourier transform so that parcevals theorem applies?

user_3426322 0
Viewed 2044

Currently, I'm considering taking an image and its spectrum. Now Parceval's theorem says that both should have equal energy. However when I try to test this on some images, this doesn't seem to be the case for the numpy real FFT function.

This is the code I'm using for my test:

import numpy as np
from PIL import Image

im = np.array(Image.open('/images/building.jpeg'))
spectral_im = np.fft.rfft2(im, axes = (0,1), norm = 'ortho')

def getNorm(im):
    return np.sum(np.abs(im))

print('Norm of the image: %d' % getNorm(im))
print('Norm of the spectrum of the image: %f' % getNorm(spectral_im))
print('Difference between norms: %f' % (getNorm(im) - getNorm(spectral_im)))

I expected that the difference between the norms would be (approximately) 0 for every image, however it differs by an order of magnitude for every image I have tried. Can anyone see what I'm doing wrong?

With the help from the answers, here is the corrected code (note the extra cast to float64, otherwise they still aren't equal):

import numpy as np
from PIL import Image

im = np.array(Image.open('/images/building.jpeg')).astype('float64')
spectral_im = np.fft.fft2(im, axes = (0,1), norm = 'ortho')

def getNorm(im):
    return np.sum(np.abs(im) ** 2)

print('Norm of the image: %d' % getNorm(im))
print('Norm of the spectrum of the image: %f' % getNorm(spectral_im))
print('Difference between norms: %f' % (getNorm(im) - getNorm(spectral_im)))
fft numpy python
2 Answers

Parceval's Theorem states that the integral over the square of the signal and the fourier transform are the same. So the getNorm function should be defined as

def getNorm(im):
    return np.sum(np.abs(im)**2)

Then there is the FFT normalization issue. You need to normalize the FFT by the image area (product of dimensions):

x = np.random.rand(321, 456)
f = np.fft.fft2(x) / np.sqrt(321 * 456)
print(np.sum(np.abs(x)**2))  # 48654.563992061871
print(np.sum(np.abs(f)**2))  # 48654.563992061878

Finally, don't use rfft to verify Parceval's theorem. The problem with rfft is that it knows the spectrum is symmetric so it skips the negative half. However, this half is missing in the integral (sum). This makes it sound like it should be off by a factor of 2, but this is not the case since the DC (mean) component is retained fully by the rfft (More details can be found here). Better use the normal FFT (fft2) and save you some trouble.

user_30051670

The standard forward FFT has an extra factor of sqrt(N) in it compared to the form used in Parceval's theorem. Taking that into account makes things work as expected:

x = np.random.rand(200, 100)
f = np.fft.fft2(x)
np.allclose(np.sum(x ** 2),
            np.sum(abs(f) ** 2 / x.size))
# True
user_29378310
Related