如何获取PSD幅值及傅里叶频谱中下一峰值的频率？

针对你提供的FFT代码，下面是获取特定峰值、PSD幅值以及下一个峰值的具体实现：

1. 提取50Hz峰值的幅值 #

通过匹配频率数组xf中最接近50Hz的位置，直接读取对应幅值：

import numpy as np
import matplotlib.pyplot as plt
import scipy.fftpack

N = 600
T = 1.0 / 800.0
x = np.linspace(0.0, N*T, N)
y = np.sin(50.0 * 2.0*np.pi*x) + 0.5*np.sin(80.0 * 2.0*np.pi*x)
yf = scipy.fftpack.fft(y)
xf = np.linspace(0.0, 1.0/(2.0*T), N//2)
fft_amp = 2.0/N * np.abs(yf[:N//2])

# 定位50Hz对应的数组索引
target_freq = 50
idx = np.argmin(np.abs(xf - target_freq))
# 获取该频率的幅值
amp_50hz = fft_amp[idx]
print(f"50Hz峰值幅值: {amp_50hz:.4f}")

# 绘图标记该点
fig, ax = plt.subplots()
ax.plot(xf, fft_amp)
ax.scatter(xf[idx], amp_50hz, color='red', label=f'50Hz ({amp_50hz:.4f})')
ax.legend()
plt.show()

2. 获取PSD（功率谱密度）的幅值 #

PSD可通过FFT结果直接计算，或用scipy.signal.welch工具实现：

方法一：基于现有FFT计算PSD #

PSD幅值等于(FFT幅值)^2 / 采样频率：

fs = 1/T  # 采样频率为800Hz
psd_amp = (fft_amp ** 2) / fs
# 提取50Hz对应的PSD幅值
psd_50hz = psd_amp[idx]
print(f"50Hz对应的PSD幅值: {psd_50hz:.6f}")

方法二：用scipy.signal.welch直接计算PSD #

from scipy import signal

f, psd = signal.welch(y, fs=fs, nperseg=N)
# 定位50Hz对应的PSD值
psd_idx = np.argmin(np.abs(f - target_freq))
print(f"Welch法计算的50Hz PSD幅值: {psd[psd_idx]:.6f}")

plt.plot(f, psd)
plt.scatter(f[psd_idx], psd[psd_idx], color='red')
plt.show()

3. 找到下一个峰值对应的频率 #

用scipy.signal.find_peaks检测所有峰值，排除50Hz后取幅值最大的那个：

from scipy.signal import find_peaks

# 检测峰值，设置最小高度过滤噪声
peaks, _ = find_peaks(fft_amp, height=0.1)
# 获取所有峰值的频率和幅值
peak_freqs = xf[peaks]
peak_amps = fft_amp[peaks]

# 按幅值从大到小排序
sorted_idx = np.argsort(peak_amps)[::-1]
sorted_peak_freqs = peak_freqs[sorted_idx]
sorted_peak_amps = peak_amps[sorted_idx]

# 第一个是50Hz，第二个即为下一个峰值
next_peak_freq = sorted_peak_freqs[1]
next_peak_amp = sorted_peak_amps[1]
print(f"下一个峰值频率: {next_peak_freq:.2f}Hz，幅值: {next_peak_amp:.4f}")

# 绘图标记两个峰值
fig, ax = plt.subplots()
ax.plot(xf, fft_amp)
ax.scatter(sorted_peak_freqs[:2], sorted_peak_amps[:2], color=['red', 'blue'], 
           label=[f'50Hz ({sorted_peak_amps[0]:.4f})', f'{next_peak_freq:.2f}Hz ({next_peak_amp:.4f})'])
ax.legend()
plt.show()

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