root/vtcross/branches/sriram/receive_path.py @ 363

#!/usr/bin/env python
#
# Copyright 2005,2006,2007,2009 Free Software Foundation, Inc.
# 
# This file is part of GNU Radio
# 
# GNU Radio is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 3, or (at your option)
# any later version.
# 
# GNU Radio is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
# 
# You should have received a copy of the GNU General Public License
# along with GNU Radio; see the file COPYING.  If not, write to
# the Free Software Foundation, Inc., 51 Franklin Street,
# Boston, MA 02110-1301, USA.
# 

from gnuradio import gr, gru, blks2
############## imports specific to spectrum sense #######################
from gnuradio import optfir, window
from gnuradio import audio
from gnuradio.eng_option import eng_option
from optparse import OptionParser
from usrpm import usrp_dbid
import math
import struct
#########################################################################
from gnuradio import usrp
from gnuradio import eng_notation
import copy
import sys
from numpy import random
# from current dir
from pick_bitrate import pick_rx_bitrate
import usrp_options
import random

# /////////////////////////////////////////////////////////////////////////////
#                              receive path
# /////////////////////////////////////////////////////////////////////////////

class receive_path(gr.hier_block2):
    def __init__(self, demod_class, rx_callback, options):

        gr.hier_block2.__init__(self, "receive_path",
                                gr.io_signature(0, 0, 0), # Input signature
                                gr.io_signature(0, 0, 0)) # Output signature

        options = copy.copy(options)    # make a copy so we can destructively modify

        self._verbose            = options.verbose
        self._rx_freq            = options.rx_freq         # receiver's center frequency
        self._rx_gain            = options.rx_gain         # receiver's gain
        self._bitrate            = options.bitrate         # desired bit rate
        self._decim              = options.decim           # Decimating rate for the USRP (prelim)
        self._samples_per_symbol = options.samples_per_symbol  # desired samples/symbol

        self._rx_callback   = rx_callback      # this callback is fired when there's a packet available
        self._demod_class   = demod_class      # the demodulator_class we're using

        if self._rx_freq is None:
            sys.stderr.write("-f FREQ or --freq FREQ or --rx-freq FREQ must be specified\n")
            raise SystemExit

        # Set up USRP source; also adjusts decim, samples_per_symbol, and bitrate
        self._setup_usrp_source(options)

        if options.show_rx_gain_range:
            print "Rx Gain Range: minimum = %g, maximum = %g, step size = %g"%tuple(self.u.gain_range())

        self.set_gain(options.rx_gain)

        # Set RF frequency
        ok = self.set_freq(self._rx_freq)
        if not ok:
            print "Failed to set Rx frequency to %s" % (eng_notation.num_to_str(self._rx_freq))
            raise ValueError, eng_notation.num_to_str(self._rx_freq)

        # copy the final answers back into options for use by demodulator
        options.samples_per_symbol = self._samples_per_symbol
        options.bitrate = self._bitrate
        options.decim = self._decim

        # Get demod_kwargs
        demod_kwargs = self._demod_class.extract_kwargs_from_options(options)

        # Design filter to get actual channel we want
        sw_decim = 1
        chan_coeffs = gr.firdes.low_pass (1.0,                  # gain
                                          sw_decim * self._samples_per_symbol, # sampling rate
                                          1.0,                  # midpoint of trans. band
                                          0.5,                  # width of trans. band
                                          gr.firdes.WIN_HANN)   # filter type 

        # Decimating channel filter
        # complex in and out, float taps
        self.chan_filt = gr.fft_filter_ccc(sw_decim, chan_coeffs)
        #self.chan_filt = gr.fir_filter_ccf(sw_decim, chan_coeffs)

        # receiver
        self.packet_receiver = \
            blks2.demod_pkts(self._demod_class(**demod_kwargs),
                             access_code=None,
                             callback=self._rx_callback,
                             threshold=-1)
    
        # Carrier Sensing Blocks
        alpha = 0.001
        thresh = 30   # in dB, will have to adjust

        if options.log_rx_power == True:
            self.probe = gr.probe_avg_mag_sqrd_cf(thresh,alpha)
            self.power_sink = gr.file_sink(gr.sizeof_float, "rxpower.dat")
            self.connect(self.chan_filt, self.probe, self.power_sink)
        else:
            self.probe = gr.probe_avg_mag_sqrd_c(thresh,alpha)
            self.connect(self.chan_filt, self.probe)

        # Display some information about the setup
        if self._verbose:
            self._print_verbage()

        #################################spectrum sense specific code#############################


        self.min_freq = 462.4825e6 # setting min and max frequency inside the init rather than taking it from the command line
        self.max_freq = 462.6425e6

        if self.min_freq > self.max_freq:
            self.min_freq, self.max_freq = self.max_freq, self.min_freq   # swap them

        self.fft_size = 512


        #if not options.real_time:
        #    realtime = False

        #else:
        #    # Attempt to enable realtime scheduling
        #    r = gr.enable_realtime_scheduling()
        #    if r == gr.RT_OK:
        #        realtime = True
        #    else:
        #        realtime = False
        #        print "Note: failed to enable realtime scheduling"

        # If the user hasn't set the fusb_* parameters on the command line,
        # pick some values that will reduce latency.

        #if 1:
        #    if options.fusb_block_size == 0 and options.fusb_nblocks == 0:
        #        if realtime:                        # be more aggressive
        #            options.fusb_block_size = gr.prefs().get_long('fusb', 'rt_block_size', 1024)
        #            options.fusb_nblocks    = gr.prefs().get_long('fusb', 'rt_nblocks', 16)
        #        else:
        #            options.fusb_block_size = gr.prefs().get_long('fusb', 'block_size', 4096)
        #            options.fusb_nblocks    = gr.prefs().get_long('fusb', 'nblocks', 16)
    
        #print "fusb_block_size =", options.fusb_block_size
        #print "fusb_nblocks    =", options.fusb_nblocks

        # build graph
        
        #self.u = usrp_options.create_usrp_source(options)


        #adc_rate = self.u.adc_rate()                # 64 MS/s
        #usrp_decim = options.decim
        #self.u.set_decim_rate(usrp_decim)
        usrp_rate = self.adc_rate / self._decim

        #self.u.set_mux(usrp.determine_rx_mux_value(self.u, options.rx_subdev_spec))
        #self.subdev = usrp.selected_subdev(self.u, options.rx_subdev_spec)
        #print "Using RX d'board %s" % (self.subdev.side_and_name(),)


        
        s2v = gr.stream_to_vector(gr.sizeof_gr_complex, self.fft_size)

        mywindow = window.blackmanharris(self.fft_size)
        fft = gr.fft_vcc(self.fft_size, True, mywindow) 
        power = 0
        for tap in mywindow:
            power += tap*tap

        
        c2mag = gr.complex_to_mag_squared(self.fft_size)
        print "print c2mag ",c2mag,"\n"
        # FIXME the log10 primitive is dog slow
        log = gr.nlog10_ff(10, self.fft_size,
                           -20*math.log10(self.fft_size)-10*math.log10(power/self.fft_size))
                
        # Set the freq_step to 75% of the actual data throughput.
        # This allows us to discard the bins on both ends of the spectrum.

        self.freq_step = 0.75 * usrp_rate
        self.min_center_freq = self.min_freq + self.freq_step/2
        nsteps = math.ceil((self.max_freq - self.min_freq) / self.freq_step)
        self.max_center_freq = self.min_center_freq + (nsteps * self.freq_step)

        self.next_freq = self.min_center_freq
        tune_delay = 1e-3
        dwell_delay = 1e-2
        tune_delay  = max(0, int(round(tune_delay * usrp_rate / self.fft_size)))  # in fft_frames
        dwell_delay = max(1, int(round(dwell_delay * usrp_rate / self.fft_size))) # in fft_frames

        self.msgq = gr.msg_queue(16)
        self._tune_callback = tune(self)        # hang on to this to keep it from being GC'd
                
        stats = gr.bin_statistics_f(self.fft_size, self.msgq,
                                    self._tune_callback, tune_delay, dwell_delay)



        # FIXME leave out the log10 until we speed it up
        #self.connect(self.u, s2v, fft, c2mag, log, stats)
        #self.connect(self.u, s2v, fft, c2mag, stats) 
        
        #if options.gain is None:
        #    # if no gain was specified, use the mid-point in dB
        #    g = self.subdev.gain_range()
        #    options.gain = float(g[0]+g[1])/2

        #self.set_gain(options.gain)
        #print "gain =", options.gain
        
        ##########################################################################################
            
        self.connect(self.u, self.chan_filt, self.packet_receiver)
        self.connect(self.u, s2v, fft, c2mag, stats) 

    def _setup_usrp_source(self, options):

        self.u = usrp_options.create_usrp_source(options)
        self.adc_rate = self.u.adc_rate()

        # derive values of bitrate, samples_per_symbol, and decim from desired info
        (self._bitrate, self._samples_per_symbol, self._decim) = \
            pick_rx_bitrate(self._bitrate, self._demod_class.bits_per_symbol(), \
                            self._samples_per_symbol, self._decim, self.adc_rate)

        self.u.set_decim(self._decim)

    def set_freq(self, target_freq):
        """
        Set the center frequency we're interested in.

        @param target_freq: frequency in Hz
        @rypte: bool

        Tuning is a two step process.  First we ask the front-end to
        tune as close to the desired frequency as it can.  Then we use
        the result of that operation and our target_frequency to
        determine the value for the digital up converter.
        """
        return self.u.set_center_freq(target_freq)

    def set_gain(self, gain):
        """
        Sets the analog gain in the USRP
        """
        return self.u.set_gain(gain)
        
    def bitrate(self):
        return self._bitrate

    def samples_per_symbol(self):
        return self._samples_per_symbol

    def decim(self):
        return self._decim

    def carrier_sensed(self):
        """
        Return True if we think carrier is present.
        """
        #return self.probe.level() > X
        return self.probe.unmuted()

    def carrier_threshold(self):
        """
        Return current setting in dB.
        """
        return self.probe.threshold()

    def set_carrier_threshold(self, threshold_in_db):
        """
        Set carrier threshold.

        @param threshold_in_db: set detection threshold
        @type threshold_in_db:  float (dB)
        """
        self.probe.set_threshold(threshold_in_db)
    
    @staticmethod
    def add_options(normal, expert):
        """
        Adds receiver-specific options to the Options Parser
        """
        add_freq_option(normal)
        if not normal.has_option("--bitrate"):
            normal.add_option("-r", "--bitrate", type="eng_float", default=None,
                              help="specify bitrate.  samples-per-symbol and interp/decim will be derived.")
        usrp_options.add_rx_options(normal, expert)
        normal.add_option("-v", "--verbose", action="store_true", default=False)
        expert.add_option("-S", "--samples-per-symbol", type="int", default=None,
                          help="set samples/symbol [default=%default]")
        expert.add_option("", "--rx-freq", type="eng_float", default=None,
                          help="set Rx frequency to FREQ [default=%default]", metavar="FREQ")
        expert.add_option("", "--log", action="store_true", default=False,
                          help="Log all parts of flow graph to files (CAUTION: lots of data)")
        expert.add_option("", "--log-rx-power", action="store_true", default=False,
                          help="Log receive signal power to file (CAUTION: lots of data)")

    def set_next_freq(self):
        target_freq = self.next_freq
        self.next_freq = self.next_freq + self.freq_step
        if self.next_freq >= self.max_center_freq:
            self.next_freq = self.min_center_freq

        if not self.set_freq(target_freq):
            print "Failed to set frequency to", target_freq

        return target_freq

    def _print_verbage(self):
        """
        Prints information about the receive path
        """
        print "\nReceive Path:"
        print "USRP %s"    % (self.u,)
        print "Rx gain:         %g"    % (self.u.gain(),)
        print "modulation:      %s"    % (self._demod_class.__name__)
        print "bitrate:         %sb/s" % (eng_notation.num_to_str(self._bitrate))
        print "samples/symbol:  %3d"   % (self._samples_per_symbol)
        print "decim:           %3d"   % (self._decim)
        print "Rx Frequency:    %s"    % (eng_notation.num_to_str(self._rx_freq))

    #Getting spectrum information
    #Arguments... threshold - Minimum power equivalent available in spectrum
#                 trials - Fft frames...should be 1 
#                 exceed count - Number of times threshold is exceeded...useless

    def get_spec_stats(self,threshold,trials):
        power_sum = 0 #sum of powers(each power value is determined by adding the 'fft square' points..these fft square points are essentially points from the PSD curve...and adding them gives us the power contained in the spectrum)
        counter = 0
        #exceed = 0
        while counter < trials:
                # Get the next message sent from the C++ code (blocking call).
                # It contains the center frequency and the mag squared of the fft
                m = parse_msg(self.msgq.delete_head())
                #print "printing mag sq of fft ",sum(m.data),"\n"
                power_sum = power_sum + sum(m.data)
                #if sum(m.data) > threshold
                #       exceed+=1
                #print sum(m.data),"\n"
                # Print center freq so we know that something is happening...
                #print m.center_freq
                counter +=1
        # FIXME do something useful with the data...
        
        
        # m.data are the mag_squared of the fft output (they are in the
        # standard order.  I.e., bin 0 == DC.)
        # You'll probably want to do the equivalent of "fftshift" on them
        # m.raw_data is a string that contains the binary floats.
        # You could write this as binary to a file.
        print "here"
        avg_power = power_sum/trials
        print "avg power ",avg_power,"\n"
        #ch = int(random.choice('17'))
        if avg_power >= threshold:
                return True
                #return 1,ch
        else:
                return False
                #return 0,"no"
        #print "printing average power ",avg_power,"\n"
        #if exceed >= exceed_count:
        #       return avg_power,True 
        #       return avg_power,True #True = primary user is present at this frequency 
        #else
        #       return avg_power,False
        #       return avg_power,False #False = primary user is NOT present at this frequency 

def add_freq_option(parser):
    """
    Hackery that has the -f / --freq option set both tx_freq and rx_freq
    """
    def freq_callback(option, opt_str, value, parser):
        parser.values.rx_freq = value
        parser.values.tx_freq = value

    if not parser.has_option('--freq'):
        parser.add_option('-f', '--freq', type="eng_float",
                          action="callback", callback=freq_callback,
                          help="set Tx and/or Rx frequency to FREQ [default=%default]",
                          metavar="FREQ")


##################### Spectrum sense specific classes##################################

class tune(gr.feval_dd):
    """
    This class allows C++ code to callback into python.
    """
    def __init__(self, tb):
        gr.feval_dd.__init__(self)
        self.tb = tb

    def eval(self, ignore):
        """
        This method is called from gr.bin_statistics_f when it wants to change
        the center frequency.  This method tunes the front end to the new center
        frequency, and returns the new frequency as its result.
        """
        try:
            # We use this try block so that if something goes wrong from here 
            # down, at least we'll have a prayer of knowing what went wrong.
            # Without this, you get a very mysterious:
            #
            #   terminate called after throwing an instance of 'Swig::DirectorMethodException'
            #   Aborted
            #
            # message on stderr.  Not exactly helpful ;)

            new_freq = self.tb.set_next_freq()
            return new_freq

        except Exception, e:
            print "tune: Exception: ", e


class parse_msg(object):
    def __init__(self, msg):
        self.center_freq = msg.arg1()
        self.vlen = int(msg.arg2())
        assert(msg.length() == self.vlen * gr.sizeof_float)

        # FIXME consider using Numarray or NumPy vector
        t = msg.to_string()
        self.raw_data = t
        self.data = struct.unpack('%df' % (self.vlen,), t)


########################################################################################