0

我可以使用以下方法从 RTL-SDR 捕获信号：

from rtlsdr import *

i = 0
signals = []
sdr = RtlSdr()

sdr.sample_rate = 2.8e6     
sdr.center_freq = 434.42e6
sdr.gain = 25

while True:
    samples = sdr.read_samples(1024*1024)
    decibel = 10*log10(var(samples))
    if decibel >= -10:
        signals.append(samples)
        i += 1
    if i == 2:
        break

如果我使用 Matplotlib 和 Seaborn 绘制信号，它们看起来像这样： 在此处输入图像描述

现在，我需要的是获得高于某个功率水平的所有峰值的坐标，例如，-20。

我发现了一个很有希望的各种Python 选项列表。但是，所有这些示例都使用了一个简单的 Numpy 数组，不同的算法可以轻松地使用它。

这是最好的尝试（因为我的猜测是我从 RTL-SDR 获得了一个复杂的信号，并且必须将其“转换”为具有实数值的数组？）：

import numpy as np
import peakutils.peak

real_sig = np.real(signals[0])
indexes = peakutils.peak.indexes(real_sig, thres=-20.0/max(real_sig), min_dist=1)
print("Peaks are: %s" % (indexes))

将这几行添加到上面的脚本中，我确实得到了一些输出，但是，首先，功率水平 -20 以上的五个峰值的值太多了。其次，这些值在给定的上下文中没有多大意义。

在此处输入图像描述

那么，我需要进行哪些更改才能获得有意义的结果，例如“Peak 1 is at 433.22 MHz”？

理想情况下，我应该得到类似的坐标Peak 1: X = 433.22, Y = -18.0，但我想一旦我知道如何获得正确的 X 值，我就可以自己弄清楚。

4

4 回答 4

1

类似于：

获取您signals的 y 值相对功率数组。

sort([x for x in signals > -20])
sort[:i]

为我拳峰。

对于频率范围：

frequency_peaks = [ frequencyspectrum[a] for a in signals[:i]]

但实际上你应该使用傅里叶变换和分箱（@hotpaw2 的回答）： 在此处输入图像描述

numpy.fft 会成功的：

https://docs.scipy.org/doc/numpy-1.13.0/reference/routines.fft.html

另见：

https://en.wikipedia.org/wiki/Nyquist%E2%80%93Shannon_sampling_theorem

为背景理论。

于 2018-03-19T19:32:40.837 回答
1

您缺少几个步骤。

您首先需要：从长度为 N 的 RTL-SDR 中选择一段复杂的 IQ 数据（您可能需要将原始 IQ 样本从无符号 8 位转换为有符号浮点），将其窗口化（von Hann 或 Hamming 等） .)，通过长度为 N 的 FFT 将其转换到频域，将 FFT 结果转换为对数幅度，并按频率标记 FFT 对数幅度结果 bin（例如数组元素），大致为

frequency(i) = sdr.center_freq + i * sdr.sample_rate / N

对于 0 到 N/2 箱

frequency(i) = sdr.center_freq - (N - i) * sdr.sample_rate / N

用于箱 N/2 到 N-1

然后，您可以沿着该对数幅度数组搜索峰值，并将频率标签应用于找到的峰值。

补充：您无法直接从 RTL-SDR 获取频域信息（如频率峰值）。你想要的峰不在那里。RTL-SDR 输出原始复数/IQ 时域样本，而不是频域数据。因此，首先您需要查找、研究并了解两者之间的区别。然后您可能会理解为什么需要 FFT（或 DFT、Goertzels 或小波等）来进行转换。

于 2018-03-19T17:53:16.987 回答
0

我相信我开始明白你们所有人都想告诉我什么...

目标仍然是重现如下图表（使用 Matplotlib 的 plt.psd() 方法创建）： 在此处输入图像描述 现在，我已经能够提出三个不同的代码段，每个都让我非常接近，但没有一个还很完美：

# Scipy version 1
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()

上面的一个产生了以下情节： 在此处输入图像描述 虽然这个图看起来还不错，但我完全不知道为什么中心峰的一部分丢失了，以及连接图两端的那条奇怪的线是从哪里来的。此外，我不知道如何显示功率电平和频率的正确值。

我的下一次尝试：

# Scipy version 2
from scipy.fftpack import fft, fftfreq
window = np.hamming(len(signal[0]))
sig_fft = fft(signal[0])
power = 20*np.log10(np.abs(sig_fft)*window)
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()

这让我得到以下结果： 在此处输入图像描述 显然，我又走上了正确的轨道，但我不知道如何在该版本的代码中应用汉明窗口（显然，我做错了）。和之前的尝试一样，我无法弄清楚如何在 x 轴上显示正确的频率。

我最后一次尝试使用 Numpy 而不是 Scipy：

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

结果是： 在此处输入图像描述 如果不是因为所有这些噪音，我会说这个可能最接近我想要的（如果没有错误应用的汉明窗，Scipy 版本 2 看起来会一样）。同样，不知道如何应用汉明窗来消除噪音。

我的问题是：

  • 如何在第二个和第三个代码段中应用汉明窗？
  • 如何在 x 轴上获得正确的频率（从 434.42-1.4 到 434.42+1.4 的值）？
  • 为什么在所有三个图中显示的信号功率水平都比原始图中（由 plt.psd() 创建）中的高得多？
于 2018-03-21T23:19:33.093 回答
0

我现在尝试了以下代码（灵感来自 hotpaw2 的答案以及我在 Google 上找到的内容）：

from numpy.fft import fft, fftshift
window = np.hamming(sdr.sample_rate/1e6+1)

plt.figure()
A = fft(window, 2048) / 25.5 # what do these values mean?
mag = np.abs(fftshift(A))
freq = np.linspace(sdr.center_freq/1e6-(sdr.sample_rate/1e6)/2, sdr.center_freq/1e6+(sdr.sample_rate/1e6)/2, len(A))
response = 20 * np.log10(mag)
#response = np.clip(response, -100, 100)
plt.plot(freq, response)
plt.title("Frequency response of Hamming window")
plt.ylabel("Magnitude [dB]")
plt.xlabel("Normalized frequency [cycles per sample]")
plt.axis("tight")
plt.show()

来源：https ://docs.scipy.org/doc/numpy-1.14.0/reference/generated/numpy.hamming.html

这导致以下（绝对无用）信号图： 在此处输入图像描述 我需要一个类似于我原始帖子中的 PSD 的图，然后我可以在其中检测峰值。

有人可以向我解释为什么我应该做所有这些 Hamming/FFT 的事情吗？

我想要的只是我的信号（由 RTL-SDR 接收）的表示，peakutils.peak.indexes() 方法接受并从中返回正确的峰值。

于 2018-03-19T18:59:40.307 回答