请给我写一个sobel算子卷积核转换到频域的代码，展示频谱和相位

Sobel算子通常用于图像处理中的边缘检测，它是由两个小的3x3卷积核构成的，分别用于检测水平和垂直方向的边缘。在频域（频率空间），我们可以利用傅里叶变换将卷积操作转化为简单的乘法运算，更直观地分析滤波的效果。以下是使用Python的NumPy库演示如何将Sobel算子转换到频域，并查看频谱和相位的一个示例：

import numpy as np
from scipy.signal import convolve2d
import matplotlib.pyplot as plt

# 定义Sobel算子
sobel_x = np.array([[1, 0, -1], [2, 0, -2], [1, 0, -1]])
sobel_y = sobel_x.T

# 将两个算子组合成二维卷积核
sobel_kernel = np.stack([sobel_x, sobel_y])

# 创建一个测试图像
image = np.zeros((50, 50))
image[25:-25, 25:-25] = 1

# 应用Sobel算子
image_sobel = convolve2d(image, sobel_kernel, mode='same')

# 计算离散傅里叶变换 (DFT)
fft_image = np.fft.fft2(image_sobel)

# 计算幅度和角度 (相位)
magnitude = np.abs(fft_image)
phase = np.angle(fft_image)

# 可视化频谱和相位
fig, axs = plt.subplots(nrows=2, ncols=2, figsize=(8, 6), sharex=True, sharey=True)
axs[0, 0].imshow(np.log(magnitude), cmap='gray', aspect='equal', origin='lower')
axs[0, 0].set_title('Magnitude Spectrum')

axs[0, 1].imshow(phase, cmap='hsv', aspect='equal', extent=(-np.pi, np.pi, -np.pi, np.pi), origin='lower')
axs[0, 1].set_title('Phase Spectrum')

axs[1, 0].imshow(magnitude, cmap='gray', aspect='auto', origin='lower')
axs[1, 0].set_title('Magnitude Image')

axs[1, 1].imshow(np.angle(fft_image), cmap='hsv', aspect='auto', extent=(-np.pi, np.pi, -np.pi, np.pi), origin='lower')
axs[1, 1].set_title('Phase Image')

for ax in axs.flat:
    ax.set(xticklabels=[], yticklabels=[])

plt.tight_layout()
plt.show()