How to get coordinates of complex RTL-SDR signal peaks?

Question

I'm able to capture signals from a RTL-SDR by using the following:

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

If I plot the signals using Matplotlib and Seaborn, they look something like this:

Now, what I would need is to get the coordinates of all peaks above a certain power level, e.g., -20.

I found a rather promising listing of the various Python options. However, all those examples use a simple Numpy array which the different algorithms can easily work with.

This was the best attempt (because my guess is that I get a complex signal from the RTL-SDR and have to "transform" it to an array with real values?):

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))

With these few lines added to the script above I do get some output, but, first, there are way too many values for the just five peaks above power level -20. And, second, the values don't make much sense in the given context.

So, what do I need to change to get meaningful results like "Peak 1 is at 433.22 MHz"?

Ideally, I should get coordinates like Peak 1: X = 433.22, Y = -18.0, but I guess I can figure that out myself once I know how to get the correct X-values.

Answer 1

You are missing several steps.

You first need to: Pick a segment of complex IQ data from the RTL-SDR of length N (you may need to convert the raw IQ samples from unsigned 8-bit to signed floating point), window it (von Hann or Hamming, etc.), convert it to the frequency domain via an FFT of length N, convert the FFT result to log magnitude, and label the FFT log magnitude result bins (e.g. array elements) by frequency, which will be roughly

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

for bins 0 to N/2

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

for bins N/2 to N-1

Then you can search along that log magnitude array for peaks, and apply the frequency label to the peaks found.

Added: You can't get frequency domain information (as in frequency peaks) directly from an RTL-SDR. The peaks you want are not there. An RTL-SDR outputs raw complex/IQ time domain samples, not frequency domain data. So first you need to look up, study, and understand the difference between the two. Then you might understand why an FFT (or DFT, or Goertzels, or wavelets, etc.) is needed to do the conversion.

Answer 2

Something similar to:

get your signals array of y value relative power.

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

for the i fist peaks.

For frequency range:

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

But really you should be using Fourrier transform and binning (@hotpaw2's answer):

numpy.fft will do the trick:

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

Also see:

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

for background theory.

Answer 3

I have now tried the following code (inspired by the answer from hotpaw2 and by what I found on 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()

Source: https://docs.scipy.org/doc/numpy-1.14.0/reference/generated/numpy.hamming.html

This results in the following (absolutely useless) signal plot: I need a plot like the PSD in my original posting, in which I can then detect the peaks.

Can someone, please, explain to my why I'm supposed to do all this Hamming/FFT stuff?

All I want is a representation of my signal (as received by the RTL-SDR) that the peakutils.peak.indexes() method accepts and returns the correct peaks from.

Answer 4

I believe I'm starting to understand what all of you were trying to tell me ...

The goal is still to reproduce a chart like the following (which was created using Matplotlib's plt.psd() method): Now, I have been able to come up with three different code segments, each of which got me pretty close, but none is perfect yet:

# 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()

The above one produces the following plot: While the figure doesn't look too bad, I have absolutely no idea why a part of the center peak is missing and where that strange line that connects both ends of the figure comes from. Also, I couldn't figure out how to display the proper values for power level and frequencies.

My next attempt:

# 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()

This gets me the following result: Clearly, I'm on the right track again but I have no idea how to apply the Hamming window in that version of the code (obviously, I did it the wrong way). And like in the previous attempt I couldn't figure out how to display the correct frequencies in the x-axis.

My last attempt uses Numpy instead of 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()

And the result is: I'd say this one probably comes closest to what I want (the Scipy version 2 would look the same without the wrongly applied hamming window), if it weren't for all that noise. Again, no idea how to apply the Hamming window to get rid of the noise.

My questions are:

• How do I apply the Hamming window in the second and third code segment?
• How do I get the proper frequencies (values from 434.42-1.4 to 434.42+1.4) on the x-axis?
• Why is the signal displayed on a significantly higher power level in all three plots than it is in the original plot (created by plt.psd())?