thesis: Add grid frequency measurement plots

lab-windows/signal_gen.ipynb

@@ -0,0 +1,298 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import struct\n",
+    "import random\n",
+    "import itertools\n",
+    "import datetime\n",
+    "import multiprocessing\n",
+    "from collections import defaultdict\n",
+    "import json\n",
+    "\n",
+    "from matplotlib import pyplot as plt\n",
+    "import matplotlib\n",
+    "import numpy as np\n",
+    "from scipy import signal as sig\n",
+    "from scipy import fftpack as fftpack\n",
+    "import ipywidgets\n",
+    "\n",
+    "import pydub\n",
+    "\n",
+    "from tqdm.notebook import tqdm\n",
+    "import colorednoise\n",
+    "\n",
+    "np.set_printoptions(linewidth=240)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib widget"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# From https://github.com/mubeta06/python/blob/master/signal_processing/sp/gold.py\n",
+    "preferred_pairs = {5:[[2],[1,2,3]], 6:[[5],[1,4,5]], 7:[[4],[4,5,6]],\n",
+    "                        8:[[1,2,3,6,7],[1,2,7]], 9:[[5],[3,5,6]], \n",
+    "                        10:[[2,5,9],[3,4,6,8,9]], 11:[[9],[3,6,9]]}\n",
+    "\n",
+    "def gen_gold(seq1, seq2):\n",
+    "    gold = [seq1, seq2]\n",
+    "    for shift in range(len(seq1)):\n",
+    "        gold.append(seq1 ^ np.roll(seq2, -shift))\n",
+    "    return gold\n",
+    "\n",
+    "def gold(n):\n",
+    "    n = int(n)\n",
+    "    if not n in preferred_pairs:\n",
+    "        raise KeyError('preferred pairs for %s bits unknown' % str(n))\n",
+    "    t0, t1 = preferred_pairs[n]\n",
+    "    (seq0, _st0), (seq1, _st1) = sig.max_len_seq(n, taps=t0), sig.max_len_seq(n, taps=t1)\n",
+    "    return gen_gold(seq0, seq1)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 49,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def modulate(data, nbits=5, pad=True):\n",
+    "    # 0, 1 -> -1, 1\n",
+    "    mask = np.array(gold(nbits))*2 - 1\n",
+    "    \n",
+    "    sel = mask[data>>1]\n",
+    "    data_lsb_centered = ((data&1)*2 - 1)\n",
+    "\n",
+    "    signal = (np.multiply(sel, np.tile(data_lsb_centered, (2**nbits-1, 1)).T).flatten() + 1) // 2\n",
+    "    if pad:\n",
+    "        return np.hstack([ np.zeros(len(mask)), signal, np.zeros(len(mask)) ])\n",
+    "    else:\n",
+    "        return signal"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 53,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def generate_noisy_signal(\n",
+    "        test_duration=32,\n",
+    "        test_nbits=5,\n",
+    "        test_decimation=10,\n",
+    "        test_signal_amplitude=200e-3,\n",
+    "        noise_level=10e-3):\n",
+    "\n",
+    "    #test_data = np.random.RandomState(seed=0).randint(0, 2 * (2**test_nbits), test_duration)\n",
+    "    #test_data = np.array([0, 1, 2, 3] * 50)\n",
+    "    test_data = np.array(range(test_duration))\n",
+    "    signal = np.repeat(modulate(test_data, test_nbits, pad=False) * 2.0 - 1, test_decimation) * test_signal_amplitude\n",
+    "    noise = colorednoise.powerlaw_psd_gaussian(1, len(signal)*2) * noise_level\n",
+    "    noise[-int(1.5*len(signal)):][:len(signal)] += signal\n",
+    "\n",
+    "    return noise+50\n",
+    "    #with open(f'mains_sim_signals/dsss_test_noisy_padded.bin', 'wb') as f:\n",
+    "    #    for e in noise:\n",
+    "    #        f.write(struct.pack('<f', e))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 54,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "application/vnd.jupyter.widget-view+json": {
+       "model_id": "c778037402024206b195a27591dc0b40",
+       "version_major": 2,
+       "version_minor": 0
+      },
+      "text/plain": [
+       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    },
+    {
+     "data": {
+      "text/plain": [
+       "[<matplotlib.lines.Line2D at 0x7f2ab8cabca0>]"
+      ]
+     },
+     "execution_count": 54,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "fig, ax = plt.subplots()\n",
+    "ax.plot(generate_noisy_signal())"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 58,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "with open('data/ref_sig_audio_test3.bin', 'wb') as f:\n",
+    "    for x in generate_noisy_signal():\n",
+    "        f.write(struct.pack('f', x))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def synthesize_sine(freqs, freqs_sampling_rate=10.0, output_sampling_rate=44100):\n",
+    "    duration = len(freqs) / freqs_sampling_rate # seconds\n",
+    "    afreq_out = np.interp(np.linspace(0, duration, int(duration*output_sampling_rate)), np.linspace(0, duration, len(freqs)), freqs)\n",
+    "    return np.sin(np.cumsum(2*np.pi * afreq_out / output_sampling_rate))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "test_sig = synthesize_sine(generate_noisy_signal())"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "application/vnd.jupyter.widget-view+json": {
+       "model_id": "5171d9dbbe5048e2b32c3cf7f7d03744",
+       "version_major": 2,
+       "version_minor": 0
+      },
+      "text/plain": [
+       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    },
+    {
+     "data": {
+      "text/plain": [
+       "[<matplotlib.lines.Line2D at 0x7f2a90476bb0>]"
+      ]
+     },
+     "execution_count": 9,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "fig, ax = plt.subplots()\n",
+    "ax.plot(test_sig[:44100])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 56,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def save_signal_flac(filename, signal, sampling_rate=44100):\n",
+    "    signal -= np.min(signal)\n",
+    "    signal /= np.max(signal)\n",
+    "    signal -= 0.5\n",
+    "    signal *= 2**16 - 1\n",
+    "    le_bytes = signal.astype(np.int16).tobytes()\n",
+    "    seg = pydub.AudioSegment(data=le_bytes, sample_width=2, frame_rate=sampling_rate, channels=1)\n",
+    "    seg.export(filename, format='flac')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 57,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "save_signal_flac('synth_sig_test_0123_02.flac', synthesize_sine(generate_noisy_signal(), freqs_sampling_rate=10.0 * 100/128, output_sampling_rate=44100))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 11,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def emulate_adc_signal(adc_bits=12, adc_offset=0.4, adc_amplitude=0.25, freq_sampling_rate=10.0, output_sampling_rate=1000, **kwargs):\n",
+    "    signal = synthesize_sine(generate_noisy_signal(), freq_sampling_rate, output_sampling_rate)\n",
+    "    signal = signal*adc_amplitude + adc_offset\n",
+    "    smin, smax = np.min(signal), np.max(signal)\n",
+    "    if smin < 0.0 or smax > 1.0:\n",
+    "        raise UserWarning('Amplitude or offset too large: Signal out of bounds with min/max [{smin}, {smax}] of ADC range')\n",
+    "    signal *= 2**adc_bits -1\n",
+    "    return signal"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 12,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def save_adc_signal(fn, signal, dtype=np.uint16):\n",
+    "    with open(fn, 'wb') as f:\n",
+    "        f.write(signal.astype(dtype).tobytes())"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 13,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "save_adc_signal('emulated_adc_readings_01.bin', emulate_adc_signal(freq_sampling_rate=10.0 * 100/128))"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "labenv",
+   "language": "python",
+   "name": "labenv"
+  },
+  "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.2"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}

ma/resources/nbexport.tplx

@@ -9,5 +9,17 @@
 
 ((*- endblock header -*))
 
-((* block maketitle *))\vspace*{3cm}((* endblock maketitle *))
+((* block maketitle *))\vspace*{1cm}((* endblock maketitle *))
+
+((* block stream *))
+    ((*- if output.name != 'stderr' -*))
+        ((( super() )))
+    ((*- endif -*))
+((* endblock stream *))
+
+((* block data_text *))
+    ((*- if 'application/vnd.jupyter.widget-view+json' not in output.data -*))
+        ((( super() )))
+    ((*- endif -*))
+((* endblock data_text *))
 

ma/safety_reset.bib

@@ -472,14 +472,14 @@
 	x-color = {#cc3300},
 	year = {2018}
 }
-
-@techreport{entsoe01,
-	author = {ENTSO-E System Protection Dynamics and WG},
-	month = mar,
-	title = {Oscillation Event 03.12.2017},
-	url = {https://docstore.entsoe.eu/Documents/SOC%20documents/Regional_Groups_Continental_Europe/OSCILLATION_REPORT_SPD.pdf},
-	year = {2018}
-}
+
+@TechReport{entsoe01,
+  author = {{ENTSO-E System Protection Dynamics and WG}},
+  title  = {Oscillation Event 03.12.2017},
+  url    = {https://docstore.entsoe.eu/Documents/SOC%20documents/Regional_Groups_Continental_Europe/OSCILLATION_REPORT_SPD.pdf},
+  month  = mar,
+  year   = {2018},
+}
 
 @article{leveson01,
 	author = {Nancy G. Leveson and Clark S. Turner},
@@ -540,14 +540,14 @@
 	title = {Smart meters in smart grid: An overview},
 	year = {2013}
 }
-
-@techreport{cenelec01,
-	author = {The CEN/CENELEC/ETSI Joint Working Group Standards Smart on for Grids},
-	month = may,
-	organization = {CEN/CENELEC/ETSI},
-	title = {Final report of the CEN/CENELEC/ETSI Joint Working Group on Standards for Smart Grids},
-	year = {2011}
-}
+
+@TechReport{cenelec01,
+  author       = {{The CEN/CENELEC/ETSI Joint Working Group Standards Smart on for Grids}},
+  title        = {Final report of the CEN/CENELEC/ETSI Joint Working Group on Standards for Smart Grids},
+  month        = may,
+  organization = {CEN/CENELEC/ETSI},
+  year         = {2011},
+}
 
 @techreport{pariente01,
 	author = {Dillon Pariente and Emmanuel Ledinot},
@@ -820,4 +820,28 @@
   year      = {2016},
 }
 
+@Article{perrin01,
+  author  = {Perrin, Trevor},
+  title   = {The Noise protocol framework, 2015},
+  journal = {URL http://noiseprotocol. org/noise. pdf},
+  year    = {2016},
+}
+
+@Article{kabalci01,
+  author = {Yasin Kabalci},
+  title  = {A survey on smart metering and smart grid communication},
+  doi    = {10.1016/j.rser.2015.12.114},
+  issn   = {1364-0321},
+  pages  = {302-318},
+  volume = {57},
+  year   = {2016},
+}
+
+@Thesis{gasior02,
+  author = {Gasior, Marek},
+  title  = {{Improving frequency resolution of discrete spectra: algorithms of three-node interpolation}},
+  url    = {https://cds.cern.ch/record/1346070},
+  year   = {2006},
+}
+
 @Comment{jabref-meta: databaseType:biblatex;}

ma/safety_reset.tex

@@ -889,7 +889,7 @@ despite numerous distortions.
 \begin{figure}
     \centering
     \includegraphics{../lab-windows/fig_out/mains_voltage_spectrum}
-    \caption{Fourier transform of an 8 hour capture of mains voltage. Data was captured using our frequency measurement
+    \caption{Fourier transform of a 24 hour capture of mains voltage. Data was captured using our frequency measurement
     sensor described in section \ref{sec-fsensor} and FFT'ed after applying a blackman window. Vertical lines indicate
     $50 \text{Hz}$ and odd harmonics.}
     \label{mains_voltage_spectrum}
@@ -1055,6 +1055,44 @@ interface and its good tolerance of system resets due to unexpected power loss.
 
 \subsection{Frequency sensor measurement results}
 
+\begin{figure}
+    \centering
+    \includegraphics{../lab-windows/fig_out/freq_meas_trace_24h}
+    \caption{Trace of grid frequency over a 24 hour window. One clearly visible feature are large positive and negative
+    transients at full hours. Times shown are UTC. Note that the european continental synchronous area that this
+    sensor is placed in covers several time zones which may result in images of daily load peaks appearing in 1 hour
+    intervals. Fig.\ \ref{freq_meas_trace_mag} contains two magnified intervals from this plot.}
+    \label{freq_meas_trace}
+\end{figure}
+\begin{figure}
+    \begin{subfigure}{\textwidth}
+        \centering
+        \includegraphics{../lab-windows/fig_out/freq_meas_trace_2h_1}
+        \caption{A 2 hour window around 00:00 UTC.}
+    \end{subfigure}
+    \begin{subfigure}{\textwidth}
+        \centering
+        \includegraphics{../lab-windows/fig_out/freq_meas_trace_2h_2}
+        \caption{A 2 hour window around 18:30 UTC.}
+    \end{subfigure}
+    \caption{Two magnified 2 hour windows of the trace from fig.\ \ref{freq_meas_trace}.}
+    \label{freq_meas_trace_mag}
+\end{figure}
+
+\begin{figure}
+    \centering
+    \includegraphics{../lab-windows/fig_out/freq_meas_spectrum}
+    \caption{Fourier transform of the 24 hour grid frequency trace in fig. \ref{freq_meas_trace} with some notable peaks
+    annotated with the corresponding period in seconds. The $\frac{1}{f}$ line indicates a pink noise spectrum. We can
+    clearly see the noise spectrum flattens below some frequency around $\frac{1}{120 \text{s}}$. This effect is due to
+    primary control actively regulating grid frequency over such time intervals. Beyond the $\frac{1}{f}$ slope starting
+    at around $1 \text{Hz}$ we can make out a white noise floor in the order of $\frac{\mu\text{Hz}}{\text{Hz}}$.
+    % TODO: where does this noise floor come from? Is it a fundamental property of the grid? Is it due to limitations of
+    % our measurement setup (such as ocxo stability/phase noise) ???
+    }
+    \label{freq_meas_spectrum}
+\end{figure}
+
 Captured raw waveform data is processed in the Jupyter Lab environment\cite{kluyver01} and grid frequency estimates are
 extracted as described in sec. \ref{frequency_estimation} using the \textcite{gasior01} technique. Appendix
 \ref{grid_freq_estimation_notebook} contains the Jupyter notebook we used for frequency measurement.
@@ -1063,6 +1101,16 @@ extracted as described in sec. \ref{frequency_estimation} using the \textcite{ga
 
 \section{Channel simulation and parameter validation}
 
+To validate all layers of our communication stack from modulation scheme to cryptography we built a prototype
+implementation in python. Implementing all components in a high-level language builds up familiartiy with the concepts
+while taking away much of the implementation complexity. For our demonstrator we will not be able to use python since
+our target platform is a cheap low-end microcontroller. Our demonstrator firmware will have to be written in a low-level
+language such as C or rust. For prototyping these languages lack flexibility compared to python.
+% FIXME introduce project outline, specs -> proto -> demo above!
+
+To validate our modulation scheme we performed a series of simulations. We produced modulated frequency data that we
+superimposed with either of simulated pink noise or an actual grid frequency measurement series.
+% FIXME do test series with simulated noise emulating measured noise spectrum
 
 \section{Implementation of a demonstrator unit} 
 
@@ -1070,6 +1118,8 @@ extracted as described in sec. \ref{frequency_estimation} using the \textcite{ga
 
 \section{Lessons learned}
 
+
+
 \chapter{Future work}
 \section{Technical standardization}
 The description of a safety reset system provided in this work could be translated into a formalized technical standard