我可以使用 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 轴上显示正确的频率值。
所以,我的问题是:
[编辑]
这是我对 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]))
如果将这些值与上图中的蓝色信号进行比较,您会发现它们非常准确。
最佳答案
我认为您的问题来自对整个信号进行 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 向你返回一个值并不意味着它是正确的或者它是你需要的值)。
由于您想要做的是在频谱中找到峰值,因此可能整个幅度的事情毕竟不是那么重要。峰值查找通常适用于频段之间的相对差异。
关于python - 如何在x轴上以正确的频率绘制信号的FFT?,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/49439510/