python用numpy库快速傅里叶变换，并绘制离散时间序列的频谱和功率谱密度图像，返回功率谱密度的三个最大值对应的频率。写出相应代码

以下是代码：

import numpy as np
import matplotlib.pyplot as plt

# 生成离散时间序列
t = np.linspace(0, 1, 1000)
x = np.sin(2 * np.pi * 10 * t) + np.sin(2 * np.pi * 20 * t) + np.sin(2 * np.pi * 30 * t)

# 进行傅里叶变换
X = np.fft.fft(x)

# 计算频谱和功率谱密度
freq = np.fft.fftfreq(len(x), t[1] - t[0])
psd = np.abs(X) ** 2 / len(x)
max_indices = np.argsort(psd)[-3:]

# 绘制频谱和功率谱密度图像
plt.figure(figsize=(8, 6))
plt.subplot(2, 1, 1)
plt.plot(freq, np.abs(X))
plt.xlabel('Frequency (Hz)')
plt.ylabel('Amplitude')
plt.title('Frequency Spectrum')
plt.subplot(2, 1, 2)
plt.plot(freq, psd)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power Spectrum Density')
plt.title('Power Spectrum Density')
plt.tight_layout()
plt.show()

# 返回功率谱密度的三个最大值对应的频率
print(freq[max_indices])

其中，`np.fft.fft(x)` 表示进行傅里叶变换，`np.abs(X)` 表示求 X 的模，`np.argsort(psd)[-3:]` 表示求 psd 的从大到小排序后的前三个索引。