Peak detection algorithm

We decided that a hueristic approach to an adaptive threshold could be using a pdf of a wider band than the one we are sensing. The idea is that we assume the noise energy is prominantly feature on the lowest part of the energy range. Therefore we needed to find the peak of the noise - mean of the main noise distribution.

After considering use of the TSpectrum class from ROOT I found a simple Matlab function that seemed to do what I wanted.

http://billauer.co.il/peakdet.html

Actually I used a Python version of the script I found posted at StackOverflow.

http://stackoverflow.com/questions/1713335/peak-finding-algorithm-for-python-scipy

Which I then converted to C++ in gnuradio. So far after playing with some of the input parameters it seems to work OK. Since we assume that the noise energy is going to be prominently featured at the low-end of the scale, this algorithm appears accurate enough to identify that noise peak.

GNU Radio energy detector implementation

After the goals of our project became more clear, I decided to revamp my energy detector implementation, particularly the 'stats' block.

The new 'stats' block is derived from gr_sync_block so we override the work() method to implement the statistics gathering. Using the old gr_bin_statistics_f as a template I arrived at a more streamlined detector based on a constant threshold value. I feel that solidifying this implementation will ease future work on an adaptive threshold.

'stats':

• for each FFT frame, calculate the average energy for a given filter bank
• accrue up to N (window size) stats
• when you have N stats gathered, generated a message with a decision
• if T > threshold, return a 1 otherwise a 0
• as new FFT frames become available, pop the oldest value and push the new one on

GNU Radio energy detector plots

Today I tested transmitting a D8PSK signal from my laptop to my desktop running the energy detector, as it is.

D8PSK

A D8PSK signal generated by benchmark_tx.py was tested first.

First, the noise + signal/noise plot:

The plot of probability of false alarm (Pfa) and probability of detection vs. threshold value:

Finally the ROC curve.

Now that the basics are in place we can look into gathering more meaningful results.

Complex Sinusoid - Lower power

A complex sinusoid signal was generated using usrp_siggen_gui.py with the following parameters:

Amplitude: 660u
TX gain: 6 dB

Complex Sinusoid - Higher power

Amplitude: 94M
TX gain: 7dB

GNU Radio signal/noise distributions

I now have my desktop setup running the detection code and my laptop using another USRP to transmit various modulated signals. After messing with the IIR alpha factor and tx-gain levels I obtained the following plots.

Again, these plots are looking at a 1MHz band centered at 895MHz. The results used 10,000 test statistics.

DQPSK:

QAM8:

GMSK:

Now the object to to calculate experimental probability of detection and false alarm.

GNU Radio the story so far...

After experimenting with a standalone C program, I decided to implement everything in GNU Radio. The hope is that eventually we can get something running in real-time, and this is the first step towards that goal.

So, at this point the 'energy detector' is implemented in two parts: ed.py and a modified version of gr_bin_statistics_f.

ed.py

This Python script follows much of the GNU Radio examples and uses the USRP as a signal source. Currently the script allows you to collect N samples, then perform an FFT followed by squaring the magnitude output. Originally the time-averaging of the PSDs was done in gr_bin_statistics_f, but the same if-not-better results were obtained using the pre-built single_pole_iir_filter. Finally the output is scaled with a log function before an output file is generated with the test statistics.

Here is the GNU Radio flow-graph:

USRP -> s2v -> FFT -> c2mag -> avg -> log -> stats

gr_bin_statistics_f

Originally, gr_bin_statistics_f was a key component in the usrp_spectrum_sense.py script. It's original function was to gather the maximum values of the FFT output and generate messages for the Python code. It implements a simple state-machine which also allows you to callback to the Python code and re-tune the USRP.

In it's current heavily-modified state, it simple calculates an average of the FFT output over a given window to generate our test stats. At this point I intend on creating a whole new block which removes much of the old functionality.

So, given a window center frequency and bandwidth, the modified code sums M bins and generates a test stat T in the form of a Python message. These messages are then tallied up by ed.py and put into an output file which is plotted directly with GNUplot.

---

In the end we have essentially the same functionality as my stand-alone C program, only now everything is implemented under GNU Radio and we are much closer to something that can be tuned to run in real-time.

libJudy: sparse dynamic arrays (associative arrays)

Seeing as libc has no hashmap component, I went searching and found Judy.

http://judy.sourceforge.net/index.html

libJudy allows you to create, among other things, associative arrays. I needed to store data with key/value pair and this fits the bill rather nicely. A set of macros have been defined to aid in the use of the library. A small summary of the basics can be found here.

Now, for associative arrays....

// first create a null pointer for the array
// which is allocated by Judy
Pvoid_t hash = (Pvoid_t) NULL;

// simple pointer to the value element
PWord_t PV;
char key[] = "hash_key";

// Judy macro to insert with a string key
JSLI(PV, hash, key);
*PV = value;

// now get the value back
JSLG(PV, hash, key);
value = *PV;

NOTES:

• Judy returns values by returning a pointer to the value, not the value itself
• The same goes to inserting values, you get a pointer to the space allocated by Judy
• Must be careful when dealing with Judy array pointers

GNU Radio PSD averaging

In search of a working energy detector, I have modified gr_bin_statistics_f to calculate a time-average of PSD frames for a given dwell delay. The intent of this is to get a cleaner signal for which we can later generate our test statistics.

Below is the 893Mhz band with the default behavior (per-bin maxima):

And here is the same band with my time-averaged PSDs (dwell = 781):

Now I need to produce a vector of test statistics to calculate our experimental probability of detection/false alarm.

Building GNU Radio on Fedora 14 (x86_64)

After installing the pre-requisites listed on the GNU Radio site (http://gnuradio.org/redmine/wiki/gnuradio/FedoraInstall) I ran into an issue while compiling:

usrp2.cc:41:33: error: type/value mismatch at argument 1 in template parameter list for ‘template class boost::weak_ptr’
usrp2.cc:41:33: error:   expected a type, got ‘usrp2::usrp2::usrp2’
usrp2.cc:43:75: error: type/value mismatch at argument 1 in template parameter list for ‘template class boost::weak_ptr’
usrp2.cc:43:75: error:   expected a type, got ‘usrp2::usrp2::usrp2’
usrp2.cc: In static member function ‘static usrp2::usrp2::sptr usrp2::usrp2::find_existing_or_make_new(const std::string&, usrp2::props*, size_t)’:
usrp2.cc:60:20: error: request for member ‘expired’ in ‘p.__gnu_cxx::__normal_iterator<_Iterator, _Container>::operator-> [with _Iterator = usrp2::usrp_table_entry*, _Container = std::vector, __gnu_cxx::__normal_iterator<_Iterator, _Container>::pointer = usrp2::usrp_table_entry*]()->usrp2::usrp_table_entry::value’, which is of non-class type ‘int’
usrp2.cc:64:31: error: no matching function for call to ‘boost::shared_ptr::shared_ptr(int&)’
/usr/include/boost/smart_ptr/shared_ptr.hpp:182:5: note: candidates are: boost::shared_ptr::shared_ptr() [with T = usrp2::usrp2]
/usr/include/boost/smart_ptr/shared_ptr.hpp:169:1: note:                 boost::shared_ptr::shared_ptr(const boost::shared_ptr&)
usrp2.cc:73:23: error: expected type-specifier
usrp2.cc:73:23: error: expected ‘)’
usrp2.cc:74:30: error: no matching function for call to ‘usrp2::usrp_table_entry::usrp_table_entry(std::string&, usrp2::usrp2::sptr&)’
usrp2.cc:43:5: note: candidates are: usrp2::usrp_table_entry::usrp_table_entry(const std::string&, int)
usrp2.cc:38:27: note:                 usrp2::usrp_table_entry::usrp_table_entry(const usrp2::usrp_table_entry&)
In file included from /usr/include/boost/shared_ptr.hpp:17:0,
                 from /home/tja/Research/gnuradio/energy/gnuradio-3.3.0/usrp2/host/include/usrp2/usrp2.h:22,
                 from usrp2.cc:23:
/usr/include/boost/smart_ptr/shared_ptr.hpp: In constructor ‘boost::shared_ptr::shared_ptr(Y*) [with Y = int, T = usrp2::usrp2]’:
usrp2.cc:73:56:   instantiated from here
/usr/include/boost/smart_ptr/shared_ptr.hpp:187:50: error: cannot convert ‘int*’ to ‘usrp2::usrp2*’ in initialization

Eric Blossom reccommended using the latest GIT branch, however there was a patch to fix this issue using the 3.3.0 tarball:

diff --git a/usrp2/host/lib/usrp2.cc b/usrp2/host/lib/usrp2.cc
index f0ee564..0842482 100644
--- a/usrp2/host/lib/usrp2.cc
+++ b/usrp2/host/lib/usrp2.cc
@@ -38,9 +38,9 @@ namespace usrp2 {
   struct usrp_table_entry {
     // inteface + normalized mac addr ("eth0:01:23:45:67:89:ab")
     std::string key;
-    boost::weak_ptr  value;
+    boost::weak_ptr  value;
 
-    usrp_table_entry(const std::string &_key, boost::weak_ptr _value)
+    usrp_table_entry(const std::string &_key, boost::weak_ptr _value)
       : key(_key), value(_value) {}
   };
 
@@ -70,7 +70,7 @@ namespace usrp2 {
     // We don't have the USRP2 we're looking for
 
     // create a new one and stick it in the table.
-    usrp2::sptr r(new usrp2::usrp2(ifc, pr, rx_bufsize));
+    usrp2::sptr r(new usrp2(ifc, pr, rx_bufsize));
     usrp_table_entry t(key, r);
     s_table.push_back(t);

After copying the patch to a file it can be applied with the following command:

patch -p1 

Python Function Tracing

Found this after a bit of googling:

def traceit(frame, event, arg):
if event == "line":
    lineno = frame.f_lineno
    fn = frame.f_globals["__file__"]
    if fn.find("fft") != -1:
        print "file %s line %d" % (fn, lineno)
return traceit

Defining a traceit method like the one above lets you examine where the Python interpreter is going after you make functions calls. The other piece is a one-liner in the main function:

if __name__ == '__main__':
    sys.settrace(traceit)
    main ()