目前,我正在考虑拍摄一幅图像及其光谱。现在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
)并省去一些麻烦。