Script for 2.3GHz beacon waterfall

waterfall_13cm.py

#!/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')