Python图像频谱分析中对数零值处理及黑屏问题解决方案

import numpy as np
import matplotlib.pyplot as plt

# 图像尺寸
M, N = 256, 256

# 生成坐标
n1 = np.arange(0, M)
n2 = np.arange(0, N)
n1, n2 = np.meshgrid(n1, n2)

# 示例函数定义 (与原问题代码一致)
def signal_function_1(n1, n2):
    return np.sin(2 * np.pi * n1 / M + 3 * np.pi * n2 / N) # 归一化频率，避免过高频率

def signal_function_2(n1, n2):
    return np.sin(4 * np.pi * n1 / M) + np.cos(6 * np.pi * n2 / N)

def signal_function_3(n1, n2):
    Y = np.zeros_like(n1, dtype=complex)
    Y[0, 5] = Y[0, N-5] = 1
    return np.real(np.fft.ifft2(Y))

def signal_function_4(n1, n2):
    Y = np.zeros_like(n1, dtype=complex)
    Y[5, 0] = Y[N-5, 0] = 1
    return np.real(np.fft.ifft2(Y))

def signal_function_5(n1, n2):
    Y = np.zeros_like(n1, dtype=complex)
    Y[5, 5] = Y[N-5, N-5] = 1
    return np.real(np.fft.ifft2(Y))

# 信号计算
signals = [
    signal_function_1(n1, n2),
    signal_function_2(n1, n2),
    signal_function_3(n1, n2),
    signal_function_4(n1, n2),
    signal_function_5(n1, n2)
]

# 计算DFT和幅度谱
dft_signals = [np.fft.fft2(s) for s in signals]
magnitude_spectrums = []

for dft_s in dft_signals:
    abs_dft_s = np.abs(dft_s)
    # 使用where参数，仅对非零幅度值进行对数计算
    # output参数用于指定结果存储的数组，这里可以直接修改abs_dft_s的副本
    # 或者创建一个新的数组，并用fill_value填充不满足条件的位置

    # 方法1: 使用where参数，并设置fill_value为0或一个非常小的负数
    # 创建一个与abs_dft_s形状相同，填充为0的数组
    magnitude_spectrum = np.zeros_like(abs_dft_s, dtype=float) 
    np.log10(abs_dft_s, out=magnitude_spectrum, where=abs_dft_s > 0)
    magnitude_spectrum *= 20

    # 也可以直接在结果数组上操作，但要确保初始值是合理的
    # 例如，如果希望零幅度对应最低分贝值（如-100dB），可以初始化为该值
    # magnitude_spectrum = np.full_like(abs_dft_s, -100.0, dtype=float)
    # np.log10(abs_dft_s, out=magnitude_spectrum, where=abs_dft_s > 0)
    # magnitude_spectrum *= 20

    magnitude_spectrums.append(magnitude_spectrum)

# 可视化 (为简洁起见，这里只展示前两个信号的修正代码)
plt.figure(figsize=(12, 8))
plt.subplot(221), plt.imshow(signals[0], cmap='gray'), plt.title('Image 1')
plt.subplot(222), plt.imshow(np.fft.fftshift(magnitude_spectrums[0])), plt.title('Spectrum 1 (Shifted)')
plt.subplot(223), plt.imshow(signals[1], cmap='gray'), plt.title('Image 2')
plt.subplot(224), plt.imshow(np.fft.fftshift(magnitude_spectrums[1])), plt.title('Spectrum 2 (Shifted)')
plt.show()

# 示例：展示一个之前可能出现黑屏的信号（如signal_function_3）
plt.figure(figsize=(12, 4))
plt.subplot(121), plt.imshow(signals[2], cmap='gray'), plt.title('Image 3')
plt.subplot(122), plt.imshow(np.fft.fftshift(magnitude_spectrums[2])), plt.title('Spectrum 3 (Shifted)')
plt.show()

# ... (前面的代码保持不变，直到计算 magnitude_spectrums) ...

magnitude_spectrums_epsilon = []

for dft_s in dft_signals:
    abs_dft_s = np.abs(dft_s)
    epsilon = 1e-10 # 一个非常小的正数
    # 对所有幅度值加上epsilon，避免log10(0)
    magnitude_spectrum = 20 * np.log10(abs_dft_s + epsilon)
    magnitude_spectrums_epsilon.append(magnitude_spectrum)

# 可视化 (使用epsilon方案的结果)
plt.figure(figsize=(12, 8))
plt.subplot(221), plt.imshow(signals[0], cmap='gray'), plt.title('Image 1 (Epsilon)')
plt.subplot(222), plt.imshow(np.fft.fftshift(magnitude_spectrums_epsilon[0])), plt.title('Spectrum 1 (Epsilon, Shifted)')
plt.subplot(223), plt.imshow(signals[1], cmap='gray'), plt.title('Image 2 (Epsilon)')
plt.subplot(224), plt.imshow(np.fft.fftshift(magnitude_spectrums_epsilon[1])), plt.title('Spectrum 2 (Epsilon, Shifted)')
plt.show()