目前,我正在考虑拍摄一幅图像及其光谱。现在Parceval的定理说两者应该有相等的能量然而,当我尝试在一些图像上测试这个时,numpy真正的fft函数似乎不是这样。
这是我用于测试的代码:

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)))

我预期每个图像之间的差异将(大约)0,但是对于我试过的每个图像来说,它的数量相差了一个数量级。有人知道我做错了什么吗?
在答案的帮助下,下面是正确的代码(请注意float64的额外转换,否则它们仍然不相等):
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)))

最佳答案

Parceval定理指出信号平方上的积分和fourier变换是相同的因此getNorm函数应该定义为

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

然后是FFT标准化问题您需要按图像区域(维度的乘积)规格化fft:
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

最后,不要使用rfft来验证Parceval定理rfft的问题是它知道频谱是对称的,所以跳过了负半部分然而,这一半在积分(和)中丢失了这使得它听起来应该关闭2倍,但事实并非如此,因为dc(mean)组件完全由rfft保留(可以找到更多细节here)。最好使用普通的fft(fft2)并省去一些麻烦。

09-10 05:18
查看更多