python - 如何在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 轴上显示正确的频率值。

所以，我的问题是:

我是否以错误的方式应用了汉宁窗，或者我怎样才能摆脱所有的噪音？

如何在绘图的 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() 绘制来自同一源的两个不同信号的结果，而蓝色信号是通过使用 plt.psd() 返回的数组提供简单的 plot() 方法产生的( Pxx 除以 sample_rate 和应用于结果的 log10)。

[编辑 4 的小补充]

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

plt.plot(freqs, power_lvls/1.1)

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

.....

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

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"选项来获取功率而不是功率密度。要获得正确的功率值，您还需要缩放功率以考虑窗口效应。对于 hann 窗口，您需要除以 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 向你返回一个值并不意味着它是正确的或者它是你需要的值)。

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

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

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

由于您想要做的是在频谱中找到峰值，因此可能整个幅度的事情毕竟不是那么重要。峰值查找通常适用于频段之间的相对差异。