BladeRF: Track FM Station and Compute Angle of Incidence

This is a quick little project. I wanted to track an FM broadcast station. I’ve got the frequency set to the local station 98.9 (on the FM dial). What I’m doing is tracking the carrier as it moves across the frequency band and extracting the phase offset produced between the two antennas due to the wave phase being slightly different at each antenna.

It works but it doesn’t work as well as you’d like. I guess due to geography like terrain, atmospheric effects, and near objects the path isn’t pointing perfectly at the transmitter. It also has a lot of noise and tends to move over time.

This might be attributed to the fact that my SNR is low. Which is likely due to my antenna since I’m out of band with it on reception.

You can find the brfharness module on my main page of this website. It’s just a simple module that automates a lot of the boilerplate for the BladeRF.

I take the measurement and write it out to a file. You can load that up, graph it out, smooth it down, and do what you like. It also prints the value to the screen. You should see it change relative to you rotating the board.

The one part I’m not sure about is where I do ant_dist / wl * np.pi * 2.0 / 4 because I couldn’t figure out why I have to divide by 4. I didn’t simplify the equation because I wanted to leave the original intact. If you figure that out, please, leave a comment below about why it is wrong, what is correct, or why it is correct like it is. I’m afraid there is a bug/mistake there.

import brfharness
from scipy.signal import windows
import numpy as np
import scipy.signal as signal

def main():
    sps = 1000000
    freq = 98.9e6
    tx_gain = None
    serials = [
        '8bba228d760342a6970a52b3d2556d99',
    ]

    cards, buffer_samps = brfharness.setup(serials, sps, freq, tx_gain)

    binhz = 5000.0
    samps = int(sps)
    fftsz = int(sps / binhz)
    fftcnt = int(samps // fftsz)
    samps = fftsz * fftcnt

    cards[0].set_rx_frequency(freq)

    pl = None

    g = []

    sos = signal.butter(8, 100e3, 'low', output='sos', fs=sps)

    fout = open('plot', 'w')

    while True:
        cards[0].clear_buffer_request_samples(samps)
        cards[0].get_pending_samples()
        _, sa, sb = cards[0].get_pending_samples()

        sa = signal.sosfiltfilt(sos, sa)
        sb = signal.sosfiltfilt(sos, sb)

        stft = signal.ShortTimeFFT.from_window(('kaiser', 12.), sps, fftsz, int(fftsz // 2), fft_mode='centered')

        fa = stft.stft(sa)
        fb = stft.stft(sb)

        # Precompute the phase angle.
        pa = fa * np.conjugate(fb)

        if pl is None:
            # Precompute the needed coefficient.
            freq_center = np.linspace(freq - sps * 0.5, freq + sps * 0.5, pa.shape[0])
            sol = 299792458.0
            ant_dist = 10 / 16 / 3.2808399
            wl = sol / freq_center
            zpl = ant_dist / wl * np.pi * 2.0 / 4
            pl = np.zeros(pa.shape, np.float64)
            for x in range(0, len(pl)):
                pl[x, :] = zpl[x]
        # Compute the geometric angle from the phase angle.
        pat = np.arcsin(np.angle(pa) / pl) / np.pi * 180.0

        track = np.argsort(np.abs(fa), axis=0)[-1, :]
        track2 = np.argsort(np.abs(fb), axis=0)[-1, :]

        track_theta = pat[track, np.arange(0, len(track))]

        g.append(np.nanmean(track_theta))

        while len(g) > 20:
            g.pop(0)

        fout.write('%s\n' % g[-1])
        fout.flush()

        print('g', np.mean(g), g[-1])

if __name__ == '__main__':
    main()