Merge "Addition of Signal Analyzer tool"
diff --git a/bandpass_recorder/README b/bandpass_recorder/README
index f796401..44c3072 100644
--- a/bandpass_recorder/README
+++ b/bandpass_recorder/README
@@ -1,5 +1,5 @@
 Interference Recorder Script 1.0
-=====================================================================
+===============================================================================
 
 This script takes in a filename as in input and opens a gui which allows
 the user to adjust the frequency, bandwidth, gain, and record time. Upon
@@ -9,13 +9,13 @@
 Usage:
 ./bandpass_recorder.py filename
 
----------------------------------------------------------------------
+-------------------------------------------------------------------------------
 Requirements:
 - USRP Software Defined Radio
 - GNURadio and UHD Driver
 - "bandpass_recorder_file.py" and "bandpass_recorder_gui.py"
 
----------------------------------------------------------------------
+-------------------------------------------------------------------------------
 Changing radio behaviour:
 The recorder behaviour can be changed by modifying "bandpass_recorder_file.py"
 or "bandpass_recorder_gui.py" in 2 ways:
@@ -27,3 +27,15 @@
 this is used by the script to specify the output filename. The script
 will output the data stream to recorder.out and it will then be renamed
 to the specified filename.
+
+===============================================================================
+Signal Analysis Tool
+
+"signal_analyzer.py" can be used to analyze recordings in order to extract
+parameters from the data files created by the recorder script.
+
+Use "./signal_analyzer.py -h" to view options
+
+Note that large window/display values can crash your system as they are very
+memory intensive. Window values should be less than or equal to 1e8 and
+display values should be less than or equal to 1e7.
diff --git a/bandpass_recorder/bandpass_recorder.py b/bandpass_recorder/bandpass_recorder.py
index fb5ae30..053df1d 100755
--- a/bandpass_recorder/bandpass_recorder.py
+++ b/bandpass_recorder/bandpass_recorder.py
@@ -1,5 +1,11 @@
 #!/usr/bin/python
-import ctypes
+"""Bandpass Recorder Tool.
+
+This script uses the URSP and records signals from the air. The user can
+specify the center-frequency, gain value, width of the recording and the
+duration of the recording.
+"""
+
 from distutils.version import StrictVersion
 import os
 import sys
@@ -8,67 +14,88 @@
 from gnuradio import gr
 from PyQt4 import Qt
 
-# Import modules from gnuradio-companion generated files.
 from bandpass_recorder_file import bandpass_recorder_file
 from bandpass_recorder_gui import bandpass_recorder_gui
+import options
 
+optspec = """
+bandpass_recorder.py [options...] filename
+--
+v,visual      Displays a GUI for viewing the region to be recorded.
+f,frequency=  Specify the center frequency to record from in MHz [2412]
+g,gain=       Specify the decibel gain of the recording. [50.0]
+w,width=      Specify the width of the interval of the recording in KHz. [20000]
+d,duration=   Specify the length of the recording in seconds. [5]
+"""
+
+
+def record(gain, freq, width, duration, filename):
+  """Record a specified signal and store it into a file.
+
+  Args:
+    gain: Amplification of the signal. (dB)
+    freq: Center Frequency of the recording (MHz).
+    width: Width of the bandpass filter. (KHz)
+    duration: Length of the recording. (s)
+    filename: Path to the file that the recording is stored into.
+  """
+  print 'test'
+  file_block = bandpass_recorder_file()
+  print 'recording... ({0} seconds)\n'.format(duration)
+  file_block.set_gain(gain)
+  file_block.set_freq(freq * 1e6)
+  file_block.set_width(width * 1e3)
+
+  file_block.start()
+  time.sleep(duration)
+
+  file_block.stop()
+  file_block.wait()
+
+  # Rename output file
+  try:
+    os.rename('recorder.out', filename)
+  except OSError as e:
+    print '\nOSError [{0}]: {1}'.format(e.errno, e.strerror)
+    print 'Recording stored in \'recorder.out\'\n'
+
+
+def main():
+  o = options.Options(optspec)
+  opt, flags, extra = o.parse(sys.argv[1:])
+
+  if len(extra) != 1:
+    o.fatal('Filename missing or extra arguments present')
+    sys.exit(1)
+
+  recording_filename = extra[0]
+
+  if StrictVersion(Qt.qVersion()) >= StrictVersion('4.5.0'):
+    prefs = gr.prefs().get_string('qtgui', 'style', 'raster')
+    Qt.QApplication.setGraphicsSystem(prefs)
+    qapp = Qt.QApplication(sys.argv)
+
+  if opt.visual:
+    gui_block = bandpass_recorder_gui()
+    gui_block.start()
+    gui_block.show()
+
+    # Define callback behaviour for the gui_block: starts the file_block and
+    # begins recording once the gui_block is closed.
+    def gui_callback():
+      gain = gui_block.get_gain()
+      freq = gui_block.get_freq() / 1000000.0
+      width = gui_block.get_width()
+      record_time = gui_block.get_record_time() / 1000.0
+
+      gui_block.stop()
+      gui_block.wait()
+      record(gain, freq, width, record_time, recording_filename)
+    qapp.connect(qapp, Qt.SIGNAL('aboutToQuit()'), gui_callback)
+    qapp.exec_()
+  else:
+    record(float(opt.gain), float(opt.frequency), float(opt.width),
+           float(opt.duration), recording_filename)
 
 if __name__ == '__main__':
-  # Checks if libX11 works properly.
-  if sys.platform.startswith('linux'):
-    try:
-      x11 = ctypes.cdll.LoadLibrary('libX11.so')
-      x11.XInitThreads()
-    except:
-      print 'Warning: failed to XInitThreads()'
-
-  # Check and parse arguments.
-  if len(sys.argv) != 2:
-    print 'usage: ./bandpass_recorder_final.py filename\n'
-    sys.exit(1)
-  if StrictVersion(Qt.qVersion()) >= StrictVersion('4.5.0'):
-    Qt.QApplication.setGraphicsSystem(gr.prefs().get_string('qtgui', 'style', 'raster'))
-  qapp = Qt.QApplication(sys.argv)
-
-  # Start the GUI
-  gui_block = bandpass_recorder_gui()
-  gui_block.start()
-  gui_block.show()
-
-  # Define callback behaviour for the gui_block: starts the file_block and begins
-  # recording once the gui_block is closed.
-  def begin_recording():
-    # Retrieve parameters
-    gain = gui_block.get_gain()
-    freq = gui_block.get_freq()
-    width = gui_block.get_width()
-    record_time = gui_block.get_record_time() / 1000.0
-
-    # Stop the gui
-    gui_block.stop()
-    gui_block.wait()
-
-    # Start the file recorder
-    file_block = bandpass_recorder_file()
-    print 'recording... ({0} seconds)\n'.format(record_time)
-
-    # Set parameters
-    file_block.set_gain(gain)
-    file_block.set_freq(freq)
-    file_block.set_width(width)
-
-    file_block.start()
-    # Run the recording for the specified time.
-    time.sleep(record_time)
-
-    file_block.stop()
-    file_block.wait()
-
-    # Rename output file
-    try:
-      os.rename('recorder.out', sys.argv[1])
-    except OSError as e:
-      print '\nOSError [{0}]: {1}'.format(e.errno, e.strerror)
-      print 'Recording stored in \'recorder.out\'\n'
-  qapp.connect(qapp, Qt.SIGNAL('aboutToQuit()'), begin_recording)
-  qapp.exec_()
+  main()
diff --git a/bandpass_recorder/bandpass_recorder_gui.py b/bandpass_recorder/bandpass_recorder_gui.py
index 29ffd63..6e2a57b 100755
--- a/bandpass_recorder/bandpass_recorder_gui.py
+++ b/bandpass_recorder/bandpass_recorder_gui.py
@@ -106,16 +106,16 @@
         self.qtgui_sink_x_0.set_update_time(1.0/10)
         self._qtgui_sink_x_0_win = sip.wrapinstance(self.qtgui_sink_x_0.pyqwidget(), Qt.QWidget)
         self.top_layout.addWidget(self._qtgui_sink_x_0_win)
-        
+
         self.qtgui_sink_x_0.enable_rf_freq(False)
-        
-        
-          
+
+
+
 
         ##################################################
         # Connections
         ##################################################
-        self.connect((self.uhd_usrp_source_0, 0), (self.qtgui_sink_x_0, 0))    
+        self.connect((self.uhd_usrp_source_0, 0), (self.qtgui_sink_x_0, 0))
 
     def closeEvent(self, event):
         self.settings = Qt.QSettings("GNU Radio", "bandpass_recorder_gui")
diff --git a/bandpass_recorder/signal_analyzer.py b/bandpass_recorder/signal_analyzer.py
new file mode 100755
index 0000000..a3a4ce5
--- /dev/null
+++ b/bandpass_recorder/signal_analyzer.py
@@ -0,0 +1,130 @@
+#!/usr/bin/python
+"""Signal Analysis Program.
+
+This program takes in a recording of an interference pattern and extracts the
+frequency, duty-cycle, and square of the signal so that the signal can be
+regenerated using parameters.
+"""
+import sys
+
+import matplotlib.pyplot as plt
+import numpy as np
+from scipy.signal import decimate
+from scipy.signal import hilbert
+from scipy.signal import symiirorder1
+
+import options
+
+SAMPLE_RATE = 30e6
+DOWNSAMPLE_FACTOR = 1000
+MAX_DUTY_PERIOD_USEC = 10000.0
+MIN_DUTY_PERIOD = SAMPLE_RATE / DOWNSAMPLE_FACTOR / MAX_DUTY_PERIOD_USEC
+optspec = """
+signal_analyzer [options...] filename
+--
+w,window=   Specify the number of samples to be analyzed. [100000000]
+d,display  Display the raw data.
+o,output=   Specify the output filename for the profile.
+"""
+
+
+def get_active_rms(raw_data, active_indices):
+  """Calculates the rms of the signal when its on."""
+  square_data = np.square(np.absolute(raw_data))
+  raw_active_indices = np.repeat(active_indices, DOWNSAMPLE_FACTOR)
+  active_values = square_data[raw_active_indices]
+
+  # TODO(sunnyd): There's a lot of values here that are not part of the active
+  # signal, need to find a way to parse them out in the future
+  return np.sqrt(np.average(active_values))
+
+
+def get_period(switch_indices):
+  """Calculates the frequency of the signal peaks."""
+  period_samples = []
+  for i in range(0, len(switch_indices) - 1):
+    period_samples.append(switch_indices[i + 1]  - switch_indices[i])
+  return np.median(period_samples) * DOWNSAMPLE_FACTOR / SAMPLE_RATE
+
+
+def get_switch_indices(active_indices):
+  # Locate the indices where the signal is on
+  off_switch_indices = []
+  on_switch_indices = []
+  for i in range(len(active_indices) - 1):
+    if active_indices[i + 1] - active_indices[i] > MIN_DUTY_PERIOD:
+      off_switch_indices.append(active_indices[i])
+      on_switch_indices.append(active_indices[i+1])
+  return (off_switch_indices, on_switch_indices)
+
+
+def save_profile(raw_data, off_switch_indices, on_switch_indices, filename):
+  """Extract the fourier transform of the signal."""
+  # TODO(sunnyd): Find an average of many peaks rather than choosing an
+  # arbitrary one.
+  cur_slice = raw_data[off_switch_indices[3] * DOWNSAMPLE_FACTOR :
+                       on_switch_indices[3] * DOWNSAMPLE_FACTOR]
+  fourier = np.fft.fft(cur_slice, int(SAMPLE_RATE))
+  fourier_down = decimate(fourier, DOWNSAMPLE_FACTOR, 2, ftype='fir')
+
+  # Filter out some of the noise
+  fourier_mean = np.average(np.absolute(fourier_down))
+
+  # TODO(sunnyd): 1.5 is a bit arbitrary, need a good metric for this. STD Dev
+  # might work but could filter out a lot of useful signal if the entire band
+  # is high
+  fourier_down[np.absolute(fourier_down) < fourier_mean * 1.5] = 0
+  np.save(filename, np.array(fourier_down))
+
+
+def get_duty_cycle(active_bitmap, switch_indices):
+  """Calculate the Duty cycle of the signal."""
+  duty_cycles = []
+  for i in range(len(switch_indices) - 1):
+    duty_cycles.append(np.average(active_bitmap[switch_indices[i] :
+                                                switch_indices[i + 1]]))
+  return np.median(duty_cycles)
+
+
+def main():
+  o = options.Options(optspec)
+  opt, _, extra = o.parse(sys.argv[1:])
+
+  window_size = int(opt.window)
+  if len(extra) != 1:
+    o.fatal('Missing filename or extra arguments')
+  if not opt.output:
+    o.fatal('Missing output option required')
+
+  f = open(extra[0], 'rb')
+  raw_data = np.fromfile(f, dtype=np.complex64, count=window_size)
+  downsampled_data = decimate(raw_data, DOWNSAMPLE_FACTOR, 2, ftype='fir')
+
+  hilbert_data = hilbert(np.absolute(downsampled_data))
+  envelope = symiirorder1(hilbert_data, 1, 0.5)
+  abs_envelope = abs(envelope)
+
+  if opt.display:
+    display_size = min(opt.window, 10e6)
+    if opt.window > 10e6:
+      sys.stderr.write('Warning: window size too large to display. Using '
+                       'default value of 10e6\n')
+    plt.plot(downsampled_data[0: display_size].real, color='red', hold=True)
+    plt.plot(abs_envelope)
+    plt.show()
+
+  # Calculate useful values for processing.
+  mean = np.average(abs_envelope)
+  active_bitmap = abs_envelope > mean
+  active_indices = active_bitmap.nonzero()[0]
+
+  off_switch_indices, on_switch_indices = get_switch_indices(active_indices)
+  duty = get_duty_cycle(active_bitmap, off_switch_indices)
+  period = get_period(off_switch_indices)
+  rms = get_active_rms(raw_data, active_bitmap)
+  save_profile(raw_data, off_switch_indices, on_switch_indices, opt.output)
+  print 'Duty-Cycle: %f, Period: %fs, Active RMS: %f' % (duty, period, rms)
+
+
+if __name__ == '__main__':
+  main()