In [1]:
%pylab inline
import scipy
import scipy.fftpack
import scipy.signal

import scipy.io

Fc = 144.375e6
Fs = 250e3

pylab.rcParams['figure.figsize'] = (15.0, 6.0)

Populating the interactive namespace from numpy and matplotlib


In [2]:
samps = [ 512000, 128000, 1024000, 16000 ]

figure(figsize=(15,4))
psd(x, Fs=Fs)
title("Raw input")
None

In [3]:
samps = [ 16000, 128000, 512000, 128000 ]
N = len(samps)

fs, b = scipy.signal.welch(x, fs=Fs, nperseg=1024)
fs = scipy.fftpack.fftshift(fs)
b = scipy.fftpack.fftshift(b)

bl = np.log(b)
dbl = diff(bl)

None

In [4]:
y = sorted(dbl)

[ y[i*len(y)/100] for i in [10, 25, 50, 75, 90] ]

Out[4]:
[-0.015214484966499242,
-0.0067623277127637493,
0.00031054897998572528,
0.008446068080306901,
0.016065986410019661]

In [5]:
il = 5*len(y)/100
ih = 95*len(y)/100
yp = y[il:ih]

my_floor = np.mean(yp) + np.std(yp)*3.0

hits = (np.array(dbl) > my_floor).nonzero()[0]+1
hits = hits[ (hits > len(bl)/10) & (hits < len(bl)*9/10) ]

[ f(yp) for f in [ np.mean, np.std] ]

Out[5]:
[0.00054897979808226087, 0.0095945428733462488]

In [6]:
# Do some janky RLEish work to find contiguous channels.
# There are better numpy-ish ways to do this, but this is legible.
runs = []

last_start = None

gap_allowed = 2 # Allow some fudge factor here

for (h0,h1) in zip(hits[:-1], hits[1:]):
if h1 - h0 < gap_allowed:
if None == last_start:
last_start = h0
else:
if None != last_start:
runs.append( (last_start, h0) )
last_start = None

plot(fs, bl)
for h0,h1 in runs:
if h1-h0 > 1:
axvspan(fs[h0-1], fs[h1+1], color="red", alpha=0.2)
xlim(fs[len(bl)/10], fs[9*len(bl)/10])

[ ((fs[h0]+fs[h1])/2, (fs[h1]-fs[h0]), fs[h0], fs[h1]) for (h0,h1) in runs if h1-h0 > 1 ]

Out[6]:
[(-4638.671875, 488.28125, -4882.8125, -4394.53125),
(-2197.265625, 3417.96875, -3906.25, -488.28125),
(10253.90625, 488.28125, 10009.765625, 10498.046875),
(12573.2421875, 732.421875, 12207.03125, 12939.453125),
(14892.578125, 488.28125, 14648.4375, 15136.71875)]

In [7]:
specgram(x[:Fs/2], Fs=Fs, NFFT=1024, noverlap=512)

for (h0,h1) in runs:
axhspan(-1+2.0*float(h0)/len(bl), -1+2.0*float(h1)/len(bl), color="grey", alpha=0.5)
None

In [8]:
def extract_channel(f_sample, f_center, f_bw, x):
"""Extract a channel of width f_bw at f_center, from f_sample

Returns Fs,y, the new sample rate and an array of complex samples
"""

# Cribbing from the GNURadio xlate, we use their block diagram

my_state = {} # cheating backchannel for debug

# Create a LPF for our target bandwidth
n_taps = 100 # XXX Total random shot in the dark
r_cutoff = float(f_bw)/f_sample
base_lpf = scipy.signal.remez(n_taps, [0, r_cutoff, r_cutoff*2, 0.5], [1,0])

# Shift this LPF up to our center frequency
fwT0 = 2.0*np.pi*f_center/f_sample

lpf = base_lpf * np.exp(1.0j*fwT0 * np.arange(0, n_taps))

dec_rate = int(f_sample / (2*f_bw))

x_filtered = scipy.signal.lfilter(lpf, 1.0, x)

dec_is = np.arange(0, len(x_filtered), dec_rate)
y = x_filtered[dec_is]
y *= np.exp(-1.0j*fwT0*dec_is)

my_state["n_taps"] = n_taps
my_state["r_cutoff"] = r_cutoff
my_state["lpf"] = lpf
my_state["base_lpf"] = base_lpf
my_state["fwT0"] = fwT0
my_state["dec_rate"] = dec_rate
# my_state["dec_is"] = dec_is
# my_state["x_filtered"] = x_filtered

return f_sample/dec_rate, y, my_state

In [9]:
fs_new, y, extract_state = extract_channel(Fs, 15000, 5000, x)

specgram(y[0:100000], Fs=fs_new)

fs_new

Out[9]:
10000.0