不同图像计算Power Spectrum时峰值位置始终相同的问题求助

复制适配了一段图像功率谱分析代码，但无论输入哪张图像（哪怕理论上峰值位置应该不同），计算得出的功率谱峰值始终固定在约0.18的频率。以下是原代码：

import matplotlib.pyplot as plt
import numpy as np
import matplotlib.mlab as mlab
import matplotlib.gridspec as gridspec
 
img = mpimg.imread("image.png")
npix = image.shape[0]
img.shape[0] == img.shape[1]
 

f = np.fft.fft2(img)
fshift = np.fft.fftshift(f)
mag = 20*np.log(np.abs(fshift))

n = int( math.ceil(img.shape[0] / 2.) * 2 )

a = np.fft.rfft(img,n, axis=0)

a = a.real*a.real + a.imag*a.imag
a = a.sum(axis=1)/a.shape[1]

f = np.fft.rfftfreq(n)

n = int( math.ceil(img.shape[1] / 2.) * 2 )

a = np.fft.rfft(img,n,axis=1)

a = a.real*a.real + a.imag*a.imag

a = a.sum(axis=0)/a.shape[0]

f = np.fft.rfftfreq(n)

plt.plot(f[1:],a[1:])

plt.ylabel( 'Amplitude' )
plt.xlabel( 'Resolution' )
plt.yscale( 'log' )
plt.legend()

问题诊断 #

• 模块缺失与变量笔误：代码未导入mpimg和math模块，且存在image变量名错误（应为img）
• 彩色图像未处理：若输入是RGB彩色图像，直接对3通道数组做FFT会导致计算混乱，需转为单通道灰度图
• 未去除直流分量：图像的直流分量（像素均值）会在0频率处产生极大功率，掩盖其他频率的真实峰值
• 变量覆盖：代码先后计算行、列方向的功率谱，但后续变量覆盖了前面的结果，最终只绘制了列方向的单一结果
• 轴标注错误：x轴标注为Resolution（分辨率）是错误的，应为Frequency（频率）

修正后的代码 #

import matplotlib.pyplot as plt
import numpy as np
from matplotlib import mpimg
import math

# 读取图像并转为灰度图
img = mpimg.imread("image.png")
# 处理彩色图像：转为单通道灰度
if len(img.shape) == 3:
    img = np.mean(img, axis=2)
# 验证图像为正方形（可选，按需保留）
assert img.shape[0] == img.shape[1], "图像需为正方形"
npix = img.shape[0]

# 去除直流分量：减去像素均值，避免低频干扰
img = img - np.mean(img)

# 计算二维FFT并移位（用于频谱可视化，可选）
f = np.fft.fft2(img)
fshift = np.fft.fftshift(f)
mag = 20 * np.log(np.abs(fshift))

# 计算行方向（垂直轴）的功率谱
n_rows = int(math.ceil(img.shape[0] / 2.) * 2)
fft_rows = np.fft.rfft(img, n_rows, axis=0)
power_rows = fft_rows.real**2 + fft_rows.imag**2
power_rows = power_rows.sum(axis=1) / power_rows.shape[1]
freq_rows = np.fft.rfftfreq(n_rows)

# 计算列方向（水平轴）的功率谱
n_cols = int(math.ceil(img.shape[1] / 2.) * 2)
fft_cols = np.fft.rfft(img, n_cols, axis=1)
power_cols = fft_cols.real**2 + fft_cols.imag**2
power_cols = power_cols.sum(axis=0) / power_cols.shape[0]
freq_cols = np.fft.rfftfreq(n_cols)

# 绘制两个方向的功率谱对比
plt.figure(figsize=(12, 6))
plt.subplot(121)
plt.plot(freq_rows[1:], power_rows[1:], label="垂直方向")
plt.ylabel('功率')
plt.xlabel('频率')
plt.yscale('log')
plt.legend()

plt.subplot(122)
plt.plot(freq_cols[1:], power_cols[1:], label="水平方向", color='orange')
plt.ylabel('功率')
plt.xlabel('频率')
plt.yscale('log')
plt.legend()

plt.tight_layout()
plt.show()

关键修正说明 #

1. 新增彩色图像转灰度逻辑，确保FFT输入为单通道数组
2. 去除直流分量，消除低频峰值对结果的干扰
3. 分别保留行、列方向的功率谱结果，避免变量覆盖导致的单一输出
4. 修正模块导入与变量名错误，保证代码可正常运行
5. 修正轴标注为正确的“频率”，符合功率谱的物理意义

相关图像 #

最终功率谱结果：

测试图像1：

测试图像2：

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