验证自然图像一维功率谱Python计算代码的实现正确性

你针对自然图像功率谱绘制需求，参考公开实现方案编写了计算代码，核心代码如下：

def power_spectrum(im_path,plot=True,transform=True):
    if type(im_path) == str:
        image = PIL.Image.open(im_path)
        if transform:
            image = torchvision.transforms.functional.center_crop(image,(256,256))
            image = np.array(torchvision.transforms.functional.rgb_to_grayscale(image))
        else:
            image = np.array(image)
    else:
        image = im_path
    
    npix = image.shape[0]

    fourier_image = np.fft.fftn(image)
    fourier_amplitudes = np.abs(fourier_image)**2

    kfreq = np.fft.fftfreq(npix) * npix
    kfreq2D = np.meshgrid(kfreq, kfreq)
    knrm = np.sqrt(kfreq2D[0]**2 + kfreq2D[1]**2)

    knrm = knrm.flatten()
    fourier_amplitudes = fourier_amplitudes.flatten()

    kbins = np.arange(0.5, npix//2+1, 1.)
    kvals = 0.5 * (kbins[1:] + kbins[:-1])
    Abins, _, _ = stats.binned_statistic(knrm, fourier_amplitudes,
                                         statistic = "mean",
                                         bins = kbins)
    Abins *= np.pi * (kbins[1:]**2 - kbins[:-1]**2)
    if plot:
        plot_power_spectra([image],[kvals],[Abins])
    return kvals,Abins

实现正确性判定 #

核心的径向平均功率谱计算逻辑和参考方案对齐，但存在4个会导致结果偏差的问题：

• 预处理缺失：直接使用0-255范围的整数灰度值做FFT，未做浮点转换、去均值操作，不同曝光的图像结果不具备可比性，且直流分量能量异常高。
• 直流分量泄漏：k=0位置的直流点（knrm=0）会被统计到第一个径向bin中，导致低波数段功率被异常拉高，对应你示例结果里最左侧的上翘段，不符合自然图像功率谱近似幂律分布的特征。
• 无窗函数抑制频谱泄漏：图像边界的不连续跳变会导致FFT出现频谱泄漏，使得高频段功率估计值偏高，标准自然图像功率谱计算通常会加二维汉明/汉宁窗抑制该效应。
• 归一化缺失：FFT结果未做像素数归一化，不同尺寸图像计算出的功率谱幅值无法直接对比，也不满足Parseval能量守恒定理。
  另外当前实现仅支持正方形输入图像，若传入非正方形图像，波数网格计算会出现维度不匹配错误。

修正方向 #

• 补全预处理：灰度图转float32类型后减去全局像素均值，可按需除以255归一化到[0,1]区间；可选加二维汉宁窗抑制频谱泄漏，窗函数代码参考：

  npix = image.shape[0]
  hanning_2d = np.hanning(npix)[:, None] * np.hanning(npix)[None, :]
  image = (image - image.mean()) * hanning_2d
  

• 屏蔽直流分量：在做分箱统计前，直接将knrm<0.5位置对应的傅里叶振幅置为0，避免直流点污染低波数统计结果。
• 补全归一化：分箱得到的Abins除以npix**4（FFT变换+窗函数对应的归一化系数），保证能量守恒，支持跨尺寸结果对比。
• 非方图适配：若需支持非正方形输入，分别对两个图像维度计算kfreq，生成波数网格时对应匹配两个维度的尺寸。

测试示例素材 #

• 示例输入ImageNet图像：
  
• 示例输出功率谱结果：
  

内容的提问来源于stack exchange，提问作者aaprasad