如何在 x 轴上绘制具有正确频率的信号的 FFT？

Question

我可以使用 Matplotlib 的 plt.psd() 方法绘制从 RTL-SDR 接收到的信号，结果如下图： 我试图实现的最终目标是检索高于某个功率水平的所有峰值的坐标，例如，-20。当我从时域接收信号时，我必须先将它们转换为频域，这由以下代码完成：

signal = []
sdr = RtlSdr()
sdr.sample_rate = 2.8e
sdr.center_freq = 434.42e6

samples = sdr.read_samples(1024*1024)
signal.append(samples)

from scipy.fftpack import fft, fftfreq
window = np.hanning(len(signal[0]))
sig_fft = fft(signal[0]*window)
power = 20*np.log10(np.abs(sig_fft))
sample_freq = fftfreq(signal[0].size, sdr.sample_rate/signal[0].size)
plt.figure(figsize=(9.84, 3.94))
plt.plot(sample_freq, power)
plt.xlabel("Frequency (MHz)")
plt.ylabel("Relative power (dB)")
plt.show()

idx = np.argmax(np.abs(sig_fft))
freq = sample_freq[idx]
peak_freq = abs(freq)
print(peak_freq)

此代码生成以下图： 我未能实现的是，首先，摆脱所有的噪音，只画一条细线，就像在 psd() 图中一样。其次，要在 x 轴上显示正确的频率值。

所以，我的问题是：

• 我是否以错误的方式应用 Hanning 窗口，或者我怎样才能消除所有噪音？
• 如何在绘图的 x 轴上获得正确的频率值？

[编辑]

这是我对 welch() 方法的尝试：

from scipy.signal import welch
sample_freq, power = welch(signal[0], sdr.sample_rate, window="hamming")
plt.figure(figsize=(9.84, 3.94))
plt.semilogy(sample_freq, power)
plt.xlabel("Frequency (MHz)")
plt.ylabel("Relative power (dB)")
plt.show()

结果： 这样我就无法在任一轴上获得正确的值。此外，中心峰值的一部分丢失了，我完全不明白，而且情节中有一条连接信号两端的恼人线。

[编辑 2]

根据 Francois Gosselin 的回答：以下代码产生的结果与 mpl.psd() 方法产生的结果最相似：

from scipy.signal import welch
corr = 1.5
sample_freq, power = welch(signal[0], fs=sdr.sample_rate, window="hann", nperseg=2048, scaling="spectrum")
sample_freq = fftshift(sample_freq)
power = fftshift(power)/corr
print(sum(power))
plt.figure(figsize=(9.84, 3.94))
plt.semilogy(sample_freq, power)
plt.xlabel("Frequency (MHz)")
plt.ylabel("Relative power (dB)")
plt.show()

现在，唯一剩下的就是弄清楚如何在相应的轴上获得正确的频率（以 MHz 为单位）和功率（以 dB 为单位）值...

[编辑 3]

使用 EDIT 2 中的代码，但使用以下行代替 plt.semilogy(...)，

plt.plot((sample_freq+sdr.center_freq)/1e6, np.log10(power))

我得到： 但是，没有必要在 plot() 方法中添加一些“额外的计算”来增强功率，不是吗？ welch() 方法不应该已经返回正确的功率级别了吗？

[编辑 4]

在尝试了您编写的所有内容后，我发现获取 plt.psd() 方法返回的频率和功率数组是最容易理解和在我的代码中使用的解决方案：

Pxx, freqs = plt.psd(signals[0], NFFT=2048, Fs=sdr.sample_rate/1e6, Fc=sdr.center_freq/1e6, scale_by_freq=True, color="green")

power_lvls = 10*log10(Pxx/(sdr.sample_rate/1e6))
plt.plot(freqs, power_lvls)
plt.show()

结果图： 有趣的是 plt.psd() 方法似乎使用与我从返回的 Pxx 数组计算它们后得到的功率级别略有不同的功率级别。绿色和红色信号是使用 plt.psd() 绘制来自同一源的两个不同信号的结果，而蓝色信号是通过提供简单的 plot() 方法以及从 plt.psd() 返回的数组产生的（ Pxx 除以 sample_rate 和应用于结果的 log10）。

[EDIT 4 的小补充]

我刚刚看到，将计算出的 power_lvls 数组中的值除以 1.1 大致使信号处于与 plt.psd() 绘制的相同的功率水平：

plt.plot(freqs, power_lvls/1.1)

现在，这可能是什么原因？默认情况下，plt.psd() 使用校正值为 1.5 的 Hanning 窗口，或者我认为。

.......

使用以下两行，我现在还可以检索多个峰的坐标：

indexes = peakutils.peak.indexes(np.array(power_lvls), thres=0.6/max(power_lvls), min_dist=120)
print("\nX: {}\n\nY: {}\n".format(freqs_1[indexes], np.array(power_lvls)[indexes]))

如果您将这些值与上图中的蓝色信号进行比较，您会发现它们非常准确。

Answer 1

我认为您的问题来自对整个信号进行 FFT，这会导致频率分辨率过高，从而导致您看到的噪声。 Matplotlib psd 将信号分解为较短的重叠块，计算每个块的 FFT 并取平均值。 Scipy 信号中的函数 welch 也可以做到这一点。您将获得以 0 Hz 为中心的频谱。然后，您可以通过将中心频率添加到频率向量来偏移返回的频率向量以获得原始频率范围。

使用 welch 时，返回的频率和功率向量不是按频率升序排序的。您必须在绘图之前对输出进行 fftshift。此外，您传递的采样频率必须是浮点数。确保使用 scaling="spectrum" 选项来获取功率而不是功率密度。要获得正确的功率值，您还需要缩放功率以考虑窗口效应。对于汉窗，您需要除以 1.5。这是一个工作示例：

from scipy.signal import welch
from numpy.fft import fftshift, fft
from numpy import arange, exp, pi, mean
import matplotlib.pyplot as plt

#generate a 1000 Hz complex tone singnal
sample_rate = 48000. #sample rate must be a float
t=arange(1024*1024)/sample_rate
signal=exp(1j*2000*pi*t)
#amplitude correction factor
corr=1.5

#calculate the psd with welch
sample_freq, power = welch(signal, fs=sample_rate, window="hann", nperseg=256, noverlap=128, scaling='spectrum')

#fftshift the output 
sample_freq=fftshift(sample_freq)
power=fftshift(power)/corr

#check that the power sum is right
print sum(power)

plt.figure(figsize=(9.84, 3.94))
plt.plot(sample_freq, power)
plt.xlabel("Frequency (MHz)")
plt.ylabel("Relative power (dB)")  
plt.show()

编辑

我看到三个原因可以解释为什么您没有获得与 Matplotlib PSD 函数相同的幅度。顺便说一句，如果您查看文档，Matplotlib PSD 具有三个返回参数：PSD、频率向量和线对象，因此如果您想要与通过 Matplotlib 函数获得的 PSD 相同的 PSD，您可以在那里获取数据。但我建议您进一步阅读以确保您知道自己在做什么（仅仅因为 Matplotlib 向您返回一个值并不意味着它是正确的或者它是您需要的值）。

1. 计算分贝的方式是错误的。首先，以分贝为单位的绘图与在对数轴上的绘图不同。分贝刻度是对数的，而您在对数图上绘制的数据仍然是线性的，因此显然您不会得到相同的值。分贝比例是相对的，这意味着您将您的值与参考值进行比较。计算分贝的方法取决于您处理的单位类型。在主要物理单位（伏特、帕斯卡、米/秒等）的情况下，如果我们假设单位是伏特，您会这样做：20*log10(V/Vref) 其中 Vref 是参考值。现在，如果您正在处理平方单位：能量、强度、功率密度等，您可以： 10*log10(P/Pref) 其中 P 是平方数量。这是因为线性域中的平方相当于对数域中的乘以 2。无论如何，在您的情况下， 10*log10 表格适用。至于参考值，它是任意的，并且通常是特定领域的约定。例如，声学中声压的国际参考是2e-5 Pa。由于参考是任意的，当你没有任何标准值可以比较时，你可以将它设置为1。

2. 其次，Matplotlib psd 的默认设置是将“scale_by_freq”选项设置为 true。这为您提供了 Hz 或 Mhz 的功率。另一方面，welch 中的频谱选项为您提供每频段的功率。因此，在 Matplotlib 中，功率除以频率范围（2.8 MHz），而在韦尔奇中，功率除以频段数（2048）。如果我取两者的分贝比，我得到 10*log10(2048/2.8)=28.6 dB，这似乎与你得到的差异非常接近。

3. 最后，您使用的校正因子会根据您想要达到的效果而有所不同。您首先需要一个校正因子的原因是因为对窗口函数引入的信号进行了修改。窗口有效地降低了信号的总能量。您需要乘以校正以获取正确的能量。加窗还通过将能量扩散到相邻频带来影响频谱。这样做的结果是修改了信号中峰值的高度。因此，如果您想要正确的峰高，则必须使用校正因子，而如果您想要正确的能量，则使用另一个。 hann 窗的两者之比为 1.5。通常，幅度校正（正确的峰高）用于显示目的，而能量校正用于总能量很重要的情况。我认为 Matplotlib PSD 是幅度校正的，而我给你的例子是能量校正的，所以那里可能有 1.5 个因子。

由于您要做的是在频谱中找到峰值，所以也许整个幅度的事情毕竟不是那么重要。峰值查找通常适用于频带之间的相对差异。