Python中二维信号傅里叶变换频谱可视化：解决log10(0)引发的显示问题

import numpy as np

# 假设 dft_signal_X 是一个二维DFT结果数组
# 例如：
# dft_signal_X = np.fft.fft2(signal_function_X(n1, n2))

# 修正后的对数幅值谱计算
# 初始化一个全零数组作为输出，确保不满足条件的位置为0
magnitude_spectrum_X = np.zeros_like(dft_signal_X, dtype=float)

# 仅在 np.abs(dft_signal_X) 不为零的位置计算对数
np.log10(np.abs(dft_signal_X), out=magnitude_spectrum_X, where=np.abs(dft_signal_X) != 0)

# 乘以20转换为分贝
magnitude_spectrum_X = 20 * magnitude_spectrum_X

import numpy as np

# 假设 dft_signal_X 是一个二维DFT结果数组
# 例如：
# dft_signal_X = np.fft.fft2(signal_function_X(n1, n2))

epsilon = 1e-10  # 定义一个非常小的常数，确保 log10 的输入不为零

# 修正后的对数幅值谱计算
magnitude_spectrum_X = 20 * np.log10(np.abs(dft_signal_X) + epsilon)

import numpy as np
import matplotlib.pyplot as plt

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

# 生成坐标 n1, n2
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((M, N), 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((M, N), dtype=complex)
    Y[5, 0] = Y[M-5, 0] = 1 # 离散频率点
    return np.real(np.fft.ifft2(Y))

def signal_function_5(n1, n2):
    Y = np.zeros((M, N), dtype=complex)
    Y[5, 5] = Y[M-5, N-5] = 1 # 离散频率点
    return np.real(np.fft.ifft2(Y))

# 计算信号
signals = {
    '1': signal_function_1(n1, n2),
    '2': signal_function_2(n1, n2),
    '3': signal_function_3(n1, n2),
    '4': signal_function_4(n1, n2),
    '5': signal_function_5(n1, n2)
}

# 计算DFT并处理对数谱
magnitude_spectrums = {}
epsilon = 1e-10 # 用于避免 log10(0)

for key, signal in signals.items():
    dft_signal = np.fft.fft2(signal)
    # 零频率居中，以便更好地可视化
    dft_shifted = np.fft.fftshift(dft_signal)

    # 使用添加微小常数的方法计算对数幅值谱
    magnitude_spectrums[key] = 20 * np.log10(np.abs(dft_shifted) + epsilon)

# 可视化
plt.figure(figsize=(15, 10))

# 信号 1 和 2
plt.subplot(241), plt.imshow(signals['1'], cmap='gray'), plt.title('Signal 1')
plt.subplot(242), plt.imshow(magnitude_spectrums['1'], cmap='jet', vmin=-100, vmax=0), plt.title('Spectrum 1') # 调整vmin/vmax以更好地显示动态范围

plt.subplot(243), plt.imshow(signals['2'], cmap='gray'), plt.title('Signal 2')
plt.subplot(244), plt.imshow(magnitude_spectrums['2'], cmap='jet', vmin=-100, vmax=0), plt.title('Spectrum 2')

# 信号 3 和 4
plt.subplot(245), plt.imshow(signals['3'], cmap='gray'), plt.title('Signal 3')
plt.subplot(246), plt.imshow(magnitude_spectrums['3'], cmap='jet', vmin=-100, vmax=0), plt.title('Spectrum 3')

plt.subplot(247), plt.imshow(signals['4'], cmap='gray'), plt.title('Signal 4')
plt.subplot(248), plt.imshow(magnitude_spectrums['4'], cmap='jet', vmin=-100, vmax=0), plt.title('Spectrum 4')

plt.tight_layout()
plt.show()

# 信号 5 (单独显示，因为只有5个)
plt.figure(figsize=(8, 4))
plt.subplot(121), plt.imshow(signals['5'], cmap='gray'), plt.title('Signal 5')
plt.subplot(122), plt.imshow(magnitude_spectrums['5'], cmap='jet', vmin=-100, vmax=0), plt.title('Spectrum 5')
plt.tight_layout()
plt.show()