KukurukuScanner.py 10.7 KB
Newer Older
Jenda's avatar
Jenda committed
1 2
from __future__ import print_function

Jenda's avatar
Jenda committed
3
import os
Jenda's avatar
Jenda committed
4 5 6 7 8 9 10 11 12
import sys
import time
import threading
import getopt
import util
import hashlib
import math
import numpy as np
import xlater
Jenda's avatar
Jenda committed
13
import struct
Jenda's avatar
Jenda committed
14
import subprocess
Jenda's avatar
Jenda committed
15
import bisect
Jenda's avatar
Jenda committed
16 17 18
from datetime import datetime
from libutil import Struct, safe_cast

Jenda's avatar
Jenda committed
19 20
from gnuradio.filter import firdes

Jenda's avatar
Jenda committed
21
def channelhelper():
Jenda's avatar
Jenda committed
22
  ChannelHelperT = Struct("channelhelper", "rotator decim rate file taps carry rotpos firpos cylen fd_r")
Jenda's avatar
Jenda committed
23
  return ChannelHelperT()
Jenda's avatar
Jenda committed
24 25 26

COMPLEX64 = 8

Jenda's avatar
scanner  
Jenda committed
27
class KukurukuScanner():
Jenda's avatar
Jenda committed
28

Jenda's avatar
scanner  
Jenda committed
29
  def __init__(self, l):
Jenda's avatar
Jenda committed
30
    self.l = l
Jenda's avatar
scanner  
Jenda committed
31
    self.conf = util.ConfReader()
Jenda's avatar
Jenda committed
32 33 34 35 36 37

  def croncmp(self, i, frag):
    if frag == "*":
      return True
    if frag[:2] == "*/":
      num = safe_cast(frag[2:], int, None)
Jenda's avatar
scanner  
Jenda committed
38
      if not num or num == 0:
Jenda's avatar
Jenda committed
39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59
        self.l.l("Invalid modulo number %s"%frag[2:], "CRIT")
        return False
      if i % num == 0:
        return True
    raw = safe_cast(frag, int, None)
    if raw is not None:
      if raw == i:
        return True
    return False

  def crontest(self, s):
    pieces = s.split(" ")
    if len(pieces) != 5:
      self.l.l("malformed cron string %s"%s, "CRIT")
      return False

    d = datetime.fromtimestamp(time.time())

    ret = self.croncmp(d.minute, pieces[0]) and self.croncmp(d.hour, pieces[1]) and\
          self.croncmp(d.day, pieces[2]) and self.croncmp(d.month, pieces[3]) and\
          self.croncmp(d.weekday() + 1, pieces[4])
Jenda's avatar
Jenda committed
60
    self.l.l("%s -> %s"%(s, ret), "DBG")
Jenda's avatar
Jenda committed
61 62 63 64 65 66 67 68 69
    return ret

  def find_cronjob(self, cronframes):
    for frame in cronframes:
      if self.crontest(frame.cron):
        return frame

  def getfn(self, f, rate):
    ds = datetime.fromtimestamp(time.time()).strftime('%Y-%m-%d-%H-%M-%S')
Jenda's avatar
scanner  
Jenda committed
70 71 72 73 74 75 76 77 78 79
    if rate is not None:
      return "%i-%s-%i"%(f, ds, rate)
    else:
      return "%i-%s"%(f, ds)

  def dump_spectrum(self, acc, filename):
    f = open(filename, "wb")
    for flt in acc:
      f.write("%f\n"%flt)
    f.close()
Jenda's avatar
Jenda committed
80 81 82 83 84 85 86

  def record_long(self, cronframe):
    ''' Record selected channels from a scheduled scan '''
    starttime = time.time()
    self.sdr.tune(cronframe.freq)
    self.l.l("cron record: tune %i"%cronframe.freq, "INFO")
    self.sdrflush()
Jenda's avatar
Jenda committed
87
    peaks = []
Jenda's avatar
Jenda committed
88
    for channel in cronframe.channels:
Jenda's avatar
scanner  
Jenda committed
89
      peaks.append([channel.freq-cronframe.freq, channel.bw, channel])
Jenda's avatar
Jenda committed
90

Jenda's avatar
scanner  
Jenda committed
91
    self.do_record(peaks, 1, None, cronframe)
Jenda's avatar
Jenda committed
92 93 94 95 96

  def scan(self, scanframe):
    ''' Find peaks in spectrum, record if specified in allow/blacklist '''
    starttime = time.time()
    self.sdr.tune(scanframe.freq)
Jenda's avatar
scanner  
Jenda committed
97
    self.sdr.set_gain(self.conf.messgain, scanframe.gain)
Jenda's avatar
Jenda committed
98
    self.sdrflush()
Jenda's avatar
scanner  
Jenda committed
99
    self.l.l("scan: tune %ik, gain %i"%(scanframe.freq/1000, scanframe.gain), "INFO")
Jenda's avatar
Jenda committed
100 101 102 103 104 105 106 107

    delta = self.conf.interval - time.time() % self.conf.interval
    nbytes = int(self.conf.rate * COMPLEX64 * delta)
    nbytes -= nbytes % COMPLEX64
    sbuf = self.pipefile.read(nbytes)

    acc = self.compute_spectrum(sbuf)

Jenda's avatar
scanner  
Jenda committed
108 109 110
    if self.conf.dumpspectrum == "always":
      self.dump_spectrum(acc, self.getfn(scanframe.freq, None)+".spectrum.txt")

Jenda's avatar
Jenda committed
111 112 113 114
    floor = sorted(acc)[int(scanframe.floor * self.conf.fftw)]

    peaks = self.find_peaks(acc, floor + scanframe.sql)

Jenda's avatar
Jenda committed
115
    peaks = self.filter_blacklist(peaks, scanframe.freq)
Jenda's avatar
Jenda committed
116

Jenda's avatar
scanner  
Jenda committed
117 118 119 120 121 122 123
    if len(peaks) == 0:
      histo = self.compute_histogram(sbuf[:COMPLEX64*self.conf.fftw])
      self.update_and_set_gain(scanframe, histo)
      self.sdr.set_gain(self.conf.messgain, scanframe.gain)
    else:
      if self.conf.dumpspectrum == "on_signal":
        self.dump_spectrum(acc, self.getfn(scanframe.freq, None)+".spectrum.txt")
Jenda's avatar
scanner  
Jenda committed
124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142

      # determine whether we stumbled upon any specified channel which has PIPE set
      peaks2 = []
      for peak in peaks:
        modified = False
        for channel in scanframe.channels:
          if channel.freq - scanframe.freq - channel.bw/2 < peak[0] and \
             channel.freq - scanframe.freq + channel.bw/2 > peak[0]:

            if channel.pipe:
              peaks2.append([channel.freq - scanframe.freq, channel.bw, channel.pipe])
            else:
              peaks2.append([channel.freq - scanframe.freq, channel.bw])
            modified = True

        if not modified:
          peaks2.append(peak)

      self.do_record(peaks2, self.conf.filtermargin, sbuf, scanframe)
Jenda's avatar
Jenda committed
143

Jenda's avatar
Jenda committed
144 145
  def filter_blacklist(self, peaks, center):
    ret = []
Jenda's avatar
scanner  
Jenda committed
146 147 148
    if not self.conf.blacklist:
      return peaks

Jenda's avatar
Jenda committed
149
    for peak in peaks:
Jenda's avatar
scanner  
Jenda committed
150 151 152
      f = peak[0]+center

      pos = bisect.bisect(self.conf.blacklist, (f, None))
Jenda's avatar
Jenda committed
153 154

      if pos >= len(self.conf.blacklist):
Jenda's avatar
scanner  
Jenda committed
155
        pos -= 1
Jenda's avatar
Jenda committed
156 157

      entry = self.conf.blacklist[pos]
Jenda's avatar
scanner  
Jenda committed
158 159

      if f < entry[0] or f > entry[1]:
Jenda's avatar
Jenda committed
160 161 162 163 164 165 166
        ret.append(peak)
        continue

      self.l.l("Removing blacklisted signal %i"%(peak[0]+center), "INFO")

    return ret

Jenda's avatar
scanner  
Jenda committed
167 168 169 170 171 172 173 174 175 176 177 178
  def update_and_set_gain(self, frame, histo):
    diff = 0
    # if there is some signal in the highest 10 bins, decrease gain
    if sum(histo[-10:]) > 0:
      diff = -1
    # if there is no signal in the upper half, increase gain
    if sum(histo[-50:]) == 0:
      diff = 1

    frame.gain += diff
    frame.gain = np.clip(frame.gain, self.conf.mingain, self.conf.maxgain)

Jenda's avatar
scanner  
Jenda committed
179
  def do_record(self, peaks, safetymargin, buf, frame):
Jenda's avatar
Jenda committed
180 181

    lastact = time.time()
Jenda's avatar
scanner  
Jenda committed
182 183
    center = frame.freq
    floor = frame.floor
Jenda's avatar
scanner  
Jenda committed
184
    stoptime = time.time() + frame.stick
Jenda's avatar
Jenda committed
185 186 187

    helpers = []
    for peak in peaks:
Jenda's avatar
Jenda committed
188
      ch = channelhelper()
Jenda's avatar
Jenda committed
189 190 191 192 193 194 195

      f = peak[0]
      w = peak[1]*safetymargin

      ch.decim = math.ceil(self.conf.rate/w)
      ch.rate = self.conf.rate/ch.decim
      ch.rotator = -float(f)/self.conf.rate * 2*math.pi
Jenda's avatar
Jenda committed
196 197
      ch.rotpos = np.zeros(2, dtype=np.float32)
      ch.rotpos[0] = 1 # start with unit vector
Jenda's avatar
Jenda committed
198 199
      taps = firdes.low_pass(1, self.conf.rate, w/2, w*self.conf.transition, firdes.WIN_HAMMING)
      ch.taps = struct.pack("=%if"%len(taps), *taps)
Jenda's avatar
Jenda committed
200 201 202
      ch.firpos = np.zeros(1, dtype=np.int32)
      ch.cylen = len(ch.taps)
      ch.carry = '\0' * ch.cylen
Jenda's avatar
Jenda committed
203

Jenda's avatar
scanner  
Jenda committed
204
      basename = self.getfn(f+center, ch.rate)
Jenda's avatar
scanner  
Jenda committed
205
      if len(peak) >= 3 and peak[2] is not None:
Jenda's avatar
Jenda committed
206
        (ch.fd_r, ch.file) = os.pipe()
Jenda's avatar
scanner  
Jenda committed
207
        cmdline = peak[2].replace("_FILENAME_", basename)
Jenda's avatar
scanner  
Jenda committed
208 209
        subprocess.Popen([cmdline], shell=True, stdin=ch.fd_r, bufsize=-1)
        self.l.l("Recording \"%s\" (PIPE), firlen %i"%(cmdline, len(taps)), "INFO")
Jenda's avatar
Jenda committed
210
      else:
Jenda's avatar
scanner  
Jenda committed
211 212 213 214
        fullfile = basename + ".cfile"
        ch.file = os.open(fullfile, os.O_WRONLY|os.O_CREAT)
        ch.fd_r = None
        self.l.l("Recording to file \"%s\", firlen %i"%(fullfile, len(taps)), "INFO")
Jenda's avatar
Jenda committed
215

Jenda's avatar
Jenda committed
216 217
      helpers.append(ch)

Jenda's avatar
scanner  
Jenda committed
218
    while True and len(helpers) > 0:
Jenda's avatar
Jenda committed
219

Jenda's avatar
scanner  
Jenda committed
220
      # read data from sdr if needed
Jenda's avatar
Jenda committed
221 222
      if buf is None:
        buf = self.pipefile.read(self.conf.bufsize * COMPLEX64)
Jenda's avatar
scanner  
Jenda committed
223

Jenda's avatar
scanner  
Jenda committed
224 225 226 227 228 229
        if frame.stickactivity:
          acc = self.compute_spectrum(buf)
          for peak in peaks:
            if self.check_activity(acc, peak, floor):
              lastact = time.time()
              self.l.l("%f has activity, continuing"%peak[0], "INFO")
Jenda's avatar
Jenda committed
230 231 232

      # xlate each channel
      for ch in helpers:
Jenda's avatar
Jenda committed
233
        # ta vec s carry je nesmysl, udelal bych workery a ty by to zjistovaly podle ch.file
Jenda's avatar
Jenda committed
234 235 236 237
        xlater.xdump(buf, len(buf), ch.carry, ch.cylen, ch.taps, len(ch.taps),
                     int(ch.decim), ch.rotator, ch.rotpos, ch.firpos, ch.file)
        ch.carry = buf[-ch.cylen:]

Jenda's avatar
scanner  
Jenda committed
238 239
      if ((not frame.stickactivity) and time.time() > stoptime) or \
        (frame.stickactivity and time.time() > lastact + frame.silencegap):
Jenda's avatar
Jenda committed
240 241 242
        self.l.l("Record stop", "INFO")
        break

Jenda's avatar
scanner  
Jenda committed
243 244 245 246
      histo = self.compute_histogram(buf[:COMPLEX64*self.conf.fftw])
      self.update_and_set_gain(frame, histo)
      self.sdr.set_gain(self.conf.messgain, frame.gain)

Jenda's avatar
Jenda committed
247 248 249 250
      buf = None

    for ch in helpers:
      os.close(ch.file)
Jenda's avatar
Jenda committed
251 252
      if ch.fd_r:
        os.close(ch.fd_r)
Jenda's avatar
Jenda committed
253 254 255 256 257 258

  def check_activity(self, acc, peak, q):
    floor = sorted(acc)[int(q * self.conf.fftw)]

    binhz = self.conf.rate/self.conf.fftw    

Jenda's avatar
scanner  
Jenda committed
259 260 261 262
    startbin = int(peak[0]/binhz - peak[1]/(2*binhz)) + self.conf.fftw/2
    stopbin  = int(peak[0]/binhz + peak[1]/(2*binhz)) + self.conf.fftw/2

    print(acc)
Jenda's avatar
Jenda committed
263 264 265

    for i in range(startbin, stopbin):
      if acc[i] > floor:
Jenda's avatar
scanner  
Jenda committed
266
        print("acc %i lvl %f floor %f"%(i,acc[i], floor))
Jenda's avatar
Jenda committed
267 268 269 270
        return True

    return False

Jenda's avatar
scanner  
Jenda committed
271 272 273 274 275 276 277 278 279
  def compute_histogram(self, sbuf):
    """ Compute histogram with 100 bins """
    acc = np.zeros(100)
    dt = np.dtype("=f4")
    buf = np.frombuffer(sbuf, dtype=dt)
    for i in range(0, len(buf)):
      acc[np.clip(int(buf[i]*99), 0, 99)] += 1
    return acc

Jenda's avatar
Jenda committed
280
  def compute_spectrum(self, sbuf):
Jenda's avatar
Jenda committed
281 282

    acc = np.zeros(self.conf.fftw)
Jenda's avatar
Jenda committed
283
    iters = 0
Jenda's avatar
Jenda committed
284
    dt = np.dtype("=c8")
Jenda's avatar
Jenda committed
285
    for i in range(0, len(sbuf)-self.conf.fftw*COMPLEX64, self.conf.fftw*self.conf.fftskip*COMPLEX64): # compute short-time FFTs, sum result to acc
Jenda's avatar
Jenda committed
286 287 288 289
      buf = np.frombuffer(sbuf, count=self.conf.fftw, dtype=dt, offset = i)

      buf = buf*self.window

Jenda's avatar
scanner  
Jenda committed
290 291 292
      fft = np.fft.fft(buf)
      fft = (np.real(fft) * np.real(fft) + np.imag(fft) * np.imag(fft))/self.conf.fftw
      fft = np.log10(fft)*10
Jenda's avatar
Jenda committed
293

Jenda's avatar
Jenda committed
294 295 296
      acc += fft
      iters += 1

Jenda's avatar
Jenda committed
297 298 299 300 301 302 303 304 305 306 307 308 309 310 311
    acc = np.divide(acc, iters)

    # FFT yields a list of positive frequencies and then negative frequencies.
    # We want it in the "natural" order.
    acc = acc.tolist()[len(acc)/2:] + acc.tolist()[:len(acc)/2]

    # smooth the result with a simple triangular filter
    acc2 = [acc[0]]

    for i in range(1, len(acc)-1):
      acc2.append((acc[i-1]*0.5 + acc[i] + acc[i+1]*0.5) / 2)

    acc2.append(acc[len(acc)-1])

    return acc2
Jenda's avatar
Jenda committed
312

Jenda's avatar
Jenda committed
313 314
  def find_peaks(self, acc, floor):
    first = -1
Jenda's avatar
Jenda committed
315

Jenda's avatar
Jenda committed
316 317 318 319 320 321 322 323 324 325 326 327 328 329
    binhz = self.conf.rate/self.conf.fftw

    minspan = self.conf.minw/binhz
    maxspan = self.conf.maxw/binhz

    peaks = []

    for i in range(1, len(acc)):
      if acc[i] > floor and acc[i-1] < floor:
        first = i
      if acc[i] < floor and acc[i-1] > floor:
        if (i-first) >= minspan and (i-first) <= maxspan:
          f = binhz*(((i+first)/2)-self.conf.fftw/2)
          w = binhz*(i-first)
Jenda's avatar
scanner  
Jenda committed
330
          peaks.append([f, w])
Jenda's avatar
Jenda committed
331 332 333
          self.l.l("signal at %f width %f"%(f, w), "INFO")

    return peaks
Jenda's avatar
Jenda committed
334 335 336 337

  def sdrflush(self):
    self.pipefile.read(self.conf.bufsize * COMPLEX64)

Jenda's avatar
Jenda committed
338 339 340
  def work(self, sdr, file_r):
    self.sdr = sdr
    self.pipefile = file_r
Jenda's avatar
Jenda committed
341

Jenda's avatar
Jenda committed
342
    self.window = np.hamming(self.conf.fftw)
Jenda's avatar
Jenda committed
343 344 345
    selframe = None
    reclen = None
    while(True):
Jenda's avatar
Jenda committed
346 347 348 349
      delta = self.conf.interval - time.time() % self.conf.interval
      if delta < self.conf.skip:
        time.sleep(delta)

Jenda's avatar
Jenda committed
350 351 352 353 354 355
      selframe = self.find_cronjob(self.conf.cronframes)

      if selframe:
        self.record_long(selframe)
        continue

Jenda's avatar
Jenda committed
356 357 358 359 360 361 362 363
      if len(self.conf.scanframes) > 0:
        # no cron job now -- pick some frame at random
        ctime = "%i"%(math.floor(time.time()) / self.conf.interval)
        idx = int(hashlib.sha256(self.conf.nonce + ctime).hexdigest(), 16) % len(self.conf.scanframes)
        scanframe = self.conf.scanframes[idx]
        self.scan(scanframe)
      else:
        time.sleep(delta)
Jenda's avatar
Jenda committed
364