import requests
import xmltodict
import numpy as np
from pylab import where
import os
import threading
import time
import sys
import signal
import datetime

#
# jamDetect.py
#
# Simple multi-threaded script to get FFTs from a receiver and run a crude signal threshold jam 
# detector. The code first runs a calibration routine to get the shape of the FFT for each 
# band. If a calibration file exists, this step is skipped.
#
# In steady-state, the script gets the FFT data as fast as possible (today limited by the receiver
# to around 5Hz). The calibration is removed, and then the median is calculated and removed - this
# removes any change in the noise floor.
#
# The script finds the maximum from the adjusted spectrum and tests against a threshold. If the
# threshold is exceeded, the full spectrum is appended (a single line) to a file, and some basic
# information is written to stdout. Optionally, all received FFTs can be saved to a file.
#
# The script uses 4 threads to process L1, L2, L5, and E6 - on an Alloy and Krypton, we got the
# expected ~5Hz FFT rate. For platforms that do not support all four bands, e.g. Argon, the
# script will determine a band isn't available and terminate the appropriate thread. The rest of
# the data processing will continue. On an Argon, I got closer to 4.5Hz maximum FFT update rate.
# It is unclear why the receiver is generating data at a slightly lower rate. During my test, there
# was plenty of reserve CPU (but maybe some processing peaks from high priority tasks have an 
# impact on the single core systems).
#
# Testing to a remote receiver in Singapore resulted in a mean FFT rate a little over 2Hz. However,
# this station has a ~184ms ping from the Bay Area.
#
# The script requests the FFT from the receiver in a CSV format. This is simpler to generate and
# parse compared to the XML output. However, there was a bug in this code that allowed repeat 
# FFTs to be output. This was fixed on 2023-11-07 (head and 5.6X ready for V5.64/6.24).
#
# To plot the data in Matlab (assume the spectrum is not inverted!)
# Spec = load('10.1.150.53-L1.log');
# FFT = Spec(:,5:end);
# binWidth = 50/2048;
# centerIdx = 1024;
# centerFreq = Spec(1,3)*1e6;
# timeTag = Spec(:,1);
# freq = ((1:2048) - centerIdx)*binWidth + centerFreq*1e-6;
# mesh(timeTag,freq,FFT')
# colormap('jet')
# colorbar
# ylabel('Frequency [MHz]');
# xlabel('Time [GPS Second]');
# zlabel('Normalized Power [dB]')
#
# Copyright Trimble Inc. 2023


RFBands = ['L1', 'L2', 'L5', 'E6']
powerThres = 10 # dB (x10 on a linear scale)
# When True gets the "raw" FFT, this is very noisy so the thresholds
# need adjusting
raw  = False
# Instead of saving data when it exceeds a threshold, save all data when this is True.
saveAll  = False

# The FFT is after the AGC, a very strong jammer can adjust the noise floor. To compensate
# for this, the script performs a noise floor adjustment based on the median value. When
# the following is "True" we save the median adjusted FFT. If this is "False" we save the
# FFT received from the receiver after converting to dB.
saveAdjusted = True

# From local network testing the throughput (FFT rate) is roughly
# the same (~5Hz) on the Linux end limited by the rate the receiver
# generates the FFTs. However, while gzip will reduce the amount
# of data passed over the network, it will approximately double the
# HTTPD load on an Alloy (5% vs 10%). Therefore, for local testing,
# disabling gzip is the best option. A remote receiver may be 
# different as the Internet bandwidth may become a driving factor.
#headers={'Accept-Encoding':'gzip'}
headers={'Accept-Encoding':'None'}

# How much data to access in the calibration phase - this only
# runs if we don't have a calibration file
calLoops = 200

# Used to control the Ctrl-C program termination
ThreadingActive = True

def signal_handler(signal,frame):
  global ThreadingActive
  print('Ctrl+C Detected - shutting down ...')
  ThreadingActive = False

  time.sleep(3)
  for num in range(len(RFBands)):
    threadHandle[num].join()

  sys.exit(0)

def getFFT(URL,session,showTimeTag=True):
    r = session.get(URL, auth=(user,password), headers=headers, proxies={'http':None}, timeout=5)
    r.raise_for_status()
    tokens = r.text.split(',')
  
    timeTag = tokens[0] # in milliseconds in the first field
    if(showTimeTag):
      with mutex:
        print(timeTag)

    data = np.array(tokens[1:2049]).astype(int)
    return(timeTag, data)

def fftWorker(bandIndex):
  band = RFBands[bandIndex]

  # Derived constants 
  loc = '/xml/dynamic/rfSpectrumAnalyzer.csv?rfBand=' + band + '&filterMode=NoFilter'
  cacheFile = IPAddr + '-' + band 
  logFile = IPAddr + '-' + band 
  if(raw):
    loc += '&showRaw=1' # Turn off all smoothing
    cacheFile += '-Raw.txt'
    logFile   += '-Raw.log'
  else:
    cacheFile += '.txt'
    logFile   += '.log'
  
  mainURL = 'http://' + IPAddr + loc

  
  # Do a full XML read of the FFT to get the spectrum details, e.g. center ferequency
  # and whether the spectrum is inverted 
  configURL  = 'http://' + IPAddr
  configURL += '/xml/dynamic/rfSpectrumAnalyzer.xml?rfBand=' + band + '&filterMode=NoFilter'

  # Use requests' Session to create a keep-alive connection.
  session = requests.Session()

  r = session.get(configURL, auth=(user,password), proxies='', timeout=5)
  r.raise_for_status()
  d = xmltodict.parse(r.text)

  numZeroes = 0
  for thisData in d['saData']['points']['y']:
    if(thisData == '0'):
      numZeroes +=1

  if(numZeroes == 2048):
    # All zeroes == band not supported
    print(band,': No FFT data - terminating thread')
    return

  centerFreq = float(d['saData']['center_freq_mhz'])
  FreqBase = centerFreq - 25
  if(d['saData']['inverted_spectrum'] == '1'):
    inverted = True
  else:
    inverted = False

  weekNum = int(d['saData']['time']['week'])
  lastTime = -1

  # When the following code block completes, this will contain the
  # Median of the 2048 point FFT over a calibration period. We'll assume this
  # is the "normal" FFT and look for deviations relative to this data. To speed the
  # processing, we create a cache based on band/IP. If this exists it will be read
  # and we'll go straight to the regular run part of the script.
  specMedian = np.zeros(2048).astype(float)

  if(os.path.isfile(cacheFile) == False):
    calData = np.zeros((calLoops,2048)).astype(int)

    with mutex:      
      print('Calibrating')
    loop = 0
    while(loop < calLoops):
      _, data = getFFT(mainURL,session)
      calData[loop,:] = data
      loop += 1

    # Get the median
    specMedian = np.median(calData,axis=0)
    specMedian = 10.0 * np.log10(specMedian) # x10 to get to dB
    with open(cacheFile,'w') as fid:
      for thisData in specMedian:
        fid.write('%.4f\n' % thisData)

    with mutex:
      print(band,np.min(specMedian), np.max(specMedian))
  else:
    idx = 0
    with open(cacheFile,'r') as fid:
      for line in fid:
        token = line.rstrip()
        specMedian[idx] = token
        idx += 1

  loop = 0
  with mutex:      
    print('Running')


  startTime = datetime.datetime.now()
  numFFTs = 0

  # Main "forever" processing loop
  while(ThreadingActive):
    try:
      timeTag, data = getFFT(mainURL,session,showTimeTag=False)
    except:
      print('HTTP Err:',mainURL)
      time.sleep(0.1)
      continue

    # With the latest firmware, we should not get any duplicates,
    # incase anything changes add some protection
    if(int(timeTag) == lastTime):
      continue

    if(int(timeTag) < lastTime):
      # The week rolled
      weekNum += 1

    lastTime = int(timeTag)

    thisFFT = 10.0 * np.log10(data.astype(float))
    if(saveAdjusted == False):
      rawFFT = np.copy(thisFFT)

    # Detrend based on the calibrated spectrum shape
    thisFFT -= specMedian

    # Because the FFT is after the AGC, a strong jammer can
    # suppress the noise floor. A normal spectrum should now
    # be relatively flat at 0dB. Find the overall median, which
    # we'll assume is a noise floor change, and adjust
    noiseFloor = np.median(thisFFT)
    # adjust
    thisFFT -= noiseFloor

    maxPower = np.max(thisFFT)
    
    # Either we requested saving all FFTs, or the threshold is exceeded
    if( (saveAll == True) or (maxPower > powerThres) ):
      # Bands are 50MHz with 2048 points  
      step = 50/2048

      if(saveAll == False):
        idx = where(thisFFT > powerThres)[0]
        num = len(idx) # How many points > than the threshold (calc is in dB)
      
        # Just get the first match
        idx_max = where(thisFFT == maxPower)[0][0]

        if(inverted):
          idx_max = 2047 - idx_max
      
        # Use a mutex for the printing to prevent occassional garbage due to 
        # collisons between the muliple threads
        with mutex:      
          print("%s:%s Max: %.3f %5.2fdB idx:%4d num > %ddB:%d" % 
                  (band,timeTag,
                  FreqBase + step*idx_max + step/2, maxPower, idx_max, powerThres, num)) 
    
      # Now save the FFT
      with open(logFile,'a') as fid:
        if(inverted):
          flag = 1
        else:
          flag = 0

        # if flag is 1, the consumer of this data needs to rotate the FFT data!!!
        fid.write('%s %d %.4f %d ' % (timeTag, weekNum, centerFreq, flag))
        if(saveAdjusted == False):
          # We save the "rawFFT" - this is the data from the receiver (converted to dB) with
          # no further processing (we haven't removed the calibration or median)
          for thisData in rawFFT:
            fid.write('%.4f ' % thisData)
        else:
          # Save the FFT adjusted FFT
          for thisData in thisFFT:
            fid.write('%.4f ' % thisData)
        fid.write('\n')

    numFFTs += 1
    delta = (datetime.datetime.now() - startTime).total_seconds()
    # Keep track of the effective FFT update rate (frequency) since the test started
    freq = numFFTs / delta

    # Periodically output a diagnostic that indicates the FFT rate we are processing
    # for each band.
    if((numFFTs % 100) == 0):
      with mutex:
        print("%s: Rate=%.3fHz" % (band,freq))


  print('Processing Complete')

if __name__ == '__main__':

  # Command line can be:
  #
  # IPAddess
  # IPAddress username password
  # --all IPAddress
  # --all IPAddress username password
  if(len(sys.argv) == 2):
    IPAddr   = sys.argv[1]
    user     = 'admin'
    password = 'password'
  elif( (len(sys.argv) == 3) and (sys.argv[1] == '--all') ):
    IPAddr   = sys.argv[2]
    user     = 'admin'
    password = 'password'
    saveAll  = True
  elif(len(sys.argv) == 4):
    IPAddr   = sys.argv[1]
    user     = sys.argv[2]
    password = sys.argv[3]
  elif( (len(sys.argv) == 5) and (sys.argv[1] == '--all') ):
    IPAddr   = sys.argv[2]
    user     = sys.argv[3]
    password = sys.argv[4]
    saveAll  = True
  else:
    print('Usage:\npython jamDetect.py [--all] IPAddress [assumes admin:password]')
    print('python jamDetect.py [--all] IPAddress username password')
    print('optional --all saves all FFT data and does not run the threshold test')
    sys.exit(1)
  
  # Setup a Ctrl-C handler
  signal.signal(signal.SIGINT, signal_handler)

  # Setup a mutex - we'll use this to protect stdout printing.
  mutex = threading.Lock()

  # Kick off a thread for each band
  threadHandle = []
  for num in range(len(RFBands)):
    t = threading.Thread(target=fftWorker, args=(num,))
    threadHandle.append(t)
    t.start()
    # 50ms sleep - gives the task time to start before we start the next
    # one. This simply reduces the chance of the debug getting mangled
    # on stdout (no longer an issue as we are using a mutex) and makes 
    # it a little easier to debug.
    time.sleep(0.05)