Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Script for 2.3GHz beacon waterfall
#!/usr/bin/env python3
import numpy as np
from scipy.fftpack import fft, fftshift
from scipy.signal import blackman
import matplotlib.pyplot as plt
import sys
# Parameters
N = 2**12
averaging = 10
samp_rate = 48000
freq = 10000
span = 3000
hz_per_bin = samp_rate/N
centre_bin = int(np.round(freq/hz_per_bin))
span_bin = int(np.ceil(span/hz_per_bin))
samples = np.memmap('/home/daniel/gqrx_20180708_100311_2320855000.wav', dtype=np.int16, mode='r', offset=0x28)
total_transforms = (len(samples)//2)//(N//2) - 1
lines = total_transforms//averaging
window = blackman(N)
try:
data = np.load('data.npz')
waterfall = data['waterfall']
lines_done = data['lines_done']
except FileNotFoundError:
waterfall = np.zeros((lines, 2*span_bin), dtype=np.float32)
lines_done = 0
done = lines_done == lines
print('Loaded {}/{} lines already computed'.format(lines_done, lines))
try:
for line in range(lines_done, lines):
print('Computing line {}/{}'.format(line + 1, lines))
transforms = np.zeros((averaging, 2*span_bin), dtype=np.float32)
for transform in range(averaging):
start = (line * averaging + transform)*N
x = np.float32(samples[start:start+2*N:2])/np.iinfo(np.int16).max # (stereo WAV, but not IQ) + 1j*np.float32(samples[start + 1:start+1+2*N:2]))/np.iinfo(np.int16).max
f = fft(x * window)[centre_bin - span_bin : centre_bin + span_bin]
transforms[transform,:] = np.real(f)**2 + np.imag(f)**2
waterfall[line,:] = 10*np.log10(np.average(transforms, axis=0)) - 20*np.log10(N)
lines_done += 1
if not done:
np.savez('data.npz', waterfall=waterfall, lines_done=lines_done)
except KeyboardInterrupt:
if not done:
print('Keyboard Interrupt. Saving work')
np.savez('data.npz', waterfall=waterfall, lines_done=lines_done)
sys.exit()
print('Plotting PNG images')
chunk = 30000
q, r = divmod(waterfall.T.shape[1], chunk)
chunks = q + int(bool(r))
for i in range(chunks):
print('Plotting image {}/{}'.format(i+1, chunks))
plt.imsave('13cm_{}.png'.format(i), waterfall[i*chunk:(i+1)*chunk,::-1].T, vmin = -125, vmax=-95, cmap='viridis', origin='bottom')
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.