Broadcast FM Noodling

Now to look at some wide FM in the broadcast band. This should be the same as before, but we need to be more aggressive about filtering out the extra cruft above the pilot tone at 19kHz.

In [17]:
%pylab inline

from gnuradio import blocks
from gnuradio import gr
from gnuradio import uhd
from gnuradio import audio
from gnuradio.filter import firdes

import scipy
import scipy.signal
Populating the interactive namespace from numpy and matplotlib

In [18]:
F_station = int(104.5e6)
Fs=int(1e6)
gain=10
N = int(5e6)

Fc = F_station - 250000 # We capture to go a bit low of the center
In [19]:
# Set up our flow graph
# Cribbing off http://gnuradio.org/notebooks/GNURadioforSimulation.html
#
tb = gr.top_block()

# Set up the USRP for input     
usrp_source = uhd.usrp_source(",", uhd.stream_args(cpu_format="fc32", channels=[0]))
usrp_source.set_samp_rate(Fs)
usrp_source.set_center_freq(Fc, 0)
usrp_source.set_gain(gain, 0)

# Let's try to flush out the first bunch of samples
skip_head = blocks.skiphead(gr.sizeof_gr_complex, 1000)

# Limit ourselves to N samples
head = blocks.head(gr.sizeof_gr_complex, N)

# And a sink to dump them into
sink = blocks.vector_sink_c()


tb.connect(usrp_source, skip_head, head, sink) # Can use the handy serial connect method here
In [20]:
tb.run()
tb.stop()
usrp_source.stop()
x = np.array(sink.data())
psd(x, Fs=Fs); None
In [21]:
specgram(xp[0:Fs], Fs=Fs); None
In [22]:
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+(0.5-r_cutoff)/4, 0.5], [1,0])
    
    # Shift this LPF up to our center frequency
    fwT0 = 2.0*np.pi*(f_center-f_bw/2)/f_sample
    
    lpf = base_lpf * np.exp(1.0j*fwT0 * np.arange(0, n_taps))
    
    dec_rate = int(f_sample / 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

def decode_quad(x, gain):
    xp = x[1:] * np.conj(x[:-1])
    retval = gain * np.arctan2(np.imag(xp), np.real(xp))
    #lpf = scipy.signal.remez(30, [0.05, 0.25, 0.27, 0.5], [1,0])
    #retval = scipy.signal.lfilter(lpf, 1.0, retval)
    return retval
 
In [23]:
Fs_new, y, _ = extract_channel(Fs, 350000, 200000, x)
specgram(y[0:Fs], Fs=Fs_new)
Fs_new
Out[23]:
200000

And now, to decode the actual audio from the band-limited FM signal. We can just use decode_quad for this.

In [24]:
w = decode_quad(y, 1.0)
In [25]:
plot(w[0:1000])
Out[25]:
[<matplotlib.lines.Line2D at 0x7f2a1cde3cd0>]
In [26]:
figure(figsize=(14,6))
specgram(w[0:30000], Fs=Fs_new); None