Beacon Hunting/find_beacons_by_fourier.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import csv\n",
    "import datetime\n",
    "import os\n",
    "import statistics\n",
    "import plotly.express as px\n",
    "import plotly.graph_objects as go\n",
    "import numpy as np\n",
    "import math"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Getcha dat data.  You could theoretically use any kind of log that includes a timestamp.  \n",
    "# For this instance, I'm using an export CSV of outbound http data.\n",
    "\n",
    "df = pd.read_csv('C:\\\\path\\\\to\\\\file.csv', header=0, parse_dates=[\"@timestamp\"])\n",
    "df.dropna(inplace=True)\n",
    "df.sort_values(by='@timestamp', ascending='True')\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Start with just a simple timeseries histogram, to see count of traffic over time\n",
    "fig=px.histogram(df, x='@timestamp')\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set custom time lengths for a sampling spacing.  Try a bunch, so we can see what we get.\n",
    "# The sample spacing is in terms of seconds per sample.  Its inverse is the sample rate, \n",
    "# in samples per second.  \n",
    "\n",
    "#Express in terms of seconds.\n",
    "sampleperiods = ['1s', '10s', '20s', '30s', '60s', '120s', '300s', '600s', '1200s', '1800s', '3600s'] \n",
    "\n",
    " # The sampling rate must be at least 2* the highest frequency we're trying to find.\n",
    " # The sampling period is 1/(sample rate).\n",
    " # Therefore, the sample period must be AT MOST 1/2 the period you are trying to find.  \n",
    " # For example, to find periodic activity with a period of 1 minute (frequency=0.0167 Hz),\n",
    " # you will need a sample period of no higher than 30 sec (frequency=0.03 Hz).\n",
    " #\n",
    " # The sample periods above can find periodicities of up to 2 hours (7200 sec).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Re-index the timestamps, so we can get counts per timeslice.\n",
    "tmpdata = df['@timestamp']\n",
    "tmpdata.index = tmpdata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Apply a Fast Fourier Transform to the numeric counts per time slice, and plot it\n",
    "for period in sampleperiods:\n",
    "    # First, get the count of http events per time slice\n",
    "    countsperperiod = tmpdata.resample(period).count()\n",
    "    \n",
    "    # Apply a Real-Input Fast Fourier Transform on the counts of items per timeslice.  This gets the\n",
    "    # signal strength (amplitude) of each present frequency.  This will become the yaxis.\n",
    "    fft = np.fft.rfft(countsperperiod)  \n",
    "    \n",
    "    # turn the current sampling period into an int\n",
    "    dvalue = int(period.rstrip(\"s\")) \n",
    "    \n",
    "    # Use the Real-Input Fast Fourier Transform Frequency\n",
    "    # rfftfreq() function to get the frequencies present in the signal.\n",
    "    # This will become the xaxis.\n",
    "    #\n",
    "    # Arguments to rfftfreq():\n",
    "    # n = window length, defined as the number of items in the current sample spacing\n",
    "    # d = sample spacing, which is defined above (also equals 1/samplerate)\n",
    "    frequencies = np.fft.rfftfreq(len(countsperperiod), dvalue)\n",
    "    \n",
    "    # Now, plot the frequency (inverse of period) vs amplitude, in Hertz\n",
    "    fig = px.line(\n",
    "        x=frequencies,\n",
    "        y=(abs(fft)), #fft is a complex number; plotting its absolute value gives the amplitude\n",
    "        labels=dict(x=\"Frequency (cycle/sec)\", y=\"HTTP GET Requests\"),\n",
    "        title=\"HTTP GETs by Frequency; Sampling Period: \" + period\n",
    "    )\n",
    "    fig.update_yaxes(\n",
    "        range=(0, 60000),\n",
    "        constrain='domain'\n",
    "    )\n",
    "    fig.show()\n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# From whatever you found as the best samnple period earlier, that shows the strongest signal spikes.\n",
    "# Set this to form 'xxs', where 'xx' is the length of the best-fit periodicity\n",
    "bestperiod = '20s' "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Identify the strongest signals present in the frequency graph (spikes). \n",
    "# First, let's resample the data and get the FFT and FFTFreq for the best sampling period we found earlier.\n",
    "countsperperiod = tmpdata.resample(bestperiod).count()\n",
    "fft = abs(np.fft.rfft(countsperperiod))\n",
    "dvalue = int(bestperiod.rstrip(\"s\")) \n",
    "frequencies = np.fft.rfftfreq(len(countsperperiod), dvalue)\n",
    " \n",
    "# Get any signal spikes over CONST * stdev over the rest of the noise.  This will be the interesting stuff to look\n",
    "# at.  The amplitudes (y-values) come from the fft array found above.\n",
    "# Find the standard deviation of the remaining data, so we can use it to find the strongest signals present.  \n",
    "# Strip off the first 10% of the frequencies found, which will remove the DC component of the signal, leaving you with \n",
    "# just the actual signal spikes.\n",
    "strippedfrequencies = frequencies[frequencies > 0.1*max(frequencies)]\n",
    "strippedfft = fft[len(fft)-strippedfrequencies.size:]\n",
    "\n",
    "stdev = np.std(strippedfft)\n",
    "mean = np.mean(strippedfft)\n",
    "threshold = mean + 4*stdev  # Adjust that multiplier as needed, higher = less results\n",
    "print(\"mean: \" + str(mean) + \", stdev: \" + str(stdev) + \", threshold: \" + str(threshold))\n",
    "\n",
    "strongsignals = []\n",
    "for signal in strippedfft:\n",
    "    if (signal > threshold): \n",
    "        print(\"adding signal: \" +str(signal))\n",
    "        strongsignals.append(signal)\n",
    "        \n",
    "# Plot the frequency data after removing the DC component\n",
    "fig = px.line(\n",
    "    x=strippedfrequencies,\n",
    "    y=(abs(strippedfft)),\n",
    "    labels=dict(x=\"Frequency (cycle/sec)\", y=\"HTTP GET Requests\"),\n",
    "    title=\"HTTP GETs by Frequency With DC Removed; Sampling Period: \" + bestperiod\n",
    ")\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# For each strong signal: find the array index from the fft array\n",
    "signalindices = []\n",
    "i = 0\n",
    "while (i < len(strongsignals)):\n",
    "    matchingindex = np.where(fft == np.float64(strongsignals[i]))[0][0]\n",
    "    signalindices.append(matchingindex)\n",
    "    i += 1\n",
    "\n",
    "# Create a new array of the same size as the fft array.  Zero it out,\n",
    "# except for the indices you just found, which are the strong signals we\n",
    "# want to find the times for.\n",
    "strongsignalfrequencies = np.zeros(len(fft))\n",
    "for index in signalindices:\n",
    "    strongsignalfrequencies[index] = frequencies[index]\n",
    "    \n",
    "strongsignalamplitudes = np.zeros(len(fft))\n",
    "for index in signalindices:\n",
    "    strongsignalamplitudes[index] = fft[index]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Graph the data in the time domain, by your chosen sampling period\n",
    "fig = px.line(\n",
    "    countsperperiod,\n",
    "    labels=dict(x=\"Timestamp\", y=\"HTTP GET Requests\"),\n",
    "    title=\"HTTP GETs By Timestamp; Sampling Period: \" + bestperiod\n",
    ")\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# De-noise the data by filtering. Make an effective bandpass filter by\n",
    "# zeroing out all the frequencies except the strong ones found above.\n",
    "# Plot just the strong signal frequencies vs their amplitudes.\n",
    "fig = px.line(\n",
    "    x=frequencies,\n",
    "    y=(abs(strongsignalamplitudes)),\n",
    "    labels=dict(x=\"Frequency (cycle/sec)\", y=\"HTTP GET Requests\"),\n",
    "        title=\"Strongest Signal Frequencies; Sampling Period: \" + bestperiod\n",
    ")\n",
    "fig.update_yaxes(\n",
    "    range=(0, 1.5*max(strongsignalamplitudes)),\n",
    "    constrain='domain'\n",
    ")\n",
    "fig.update_xaxes(\n",
    "    \n",
    ")\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use the Inverse FFT to flip just the strong signals back to time-domain\n",
    "inversefft = np.fft.irfft(strongsignalamplitudes, len(countsperperiod))\n",
    "\n",
    "fig = px.line(\n",
    "    x=countsperperiod.to_frame().index,\n",
    "    y=inversefft,\n",
    "    labels=dict(x=\"Timestamp\", y=\"Count of HTTP Get Requests\"),\n",
    "        title=\"Periodic Signal\"\n",
    ")\n",
    "\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# OK.  Now, for each of our strong signals, we need to identify domains from our original data set that \n",
    "# had a count of GET requests \"near\" our signal strengths.  (It won't be spot-on, due to sample frequency \n",
    "# bin width and signal jitter.)  This will be the shortlist of domains for further investigation.\n",
    "shortlist = []\n",
    "newdf = df.groupby(['host']).size().reset_index(name='counts')\n",
    "for amplitude in strongsignals:\n",
    "    shortlist.append(newdf[ (newdf['counts'] > (amplitude*0.8)) & (newdf['counts'] < (amplitude*1.2)) ])\n",
    "    \n",
    "results = pd.concat(shortlist, ignore_index=True)\n",
    "print(results)\n",
    "results[['host','counts']].to_csv('c:\\\\hunting\\\\output.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}