From 6a49b8399fb156d427b5d15746de877ac933ba3d Mon Sep 17 00:00:00 2001 From: jaseg Date: Sun, 12 Jan 2020 20:06:55 +0100 Subject: meas ipy: some parameter tuning --- Phase Measurement Prototype.ipynb | 1803 +++++++++++++++++++++++++++++++++++++ Phase measurement Prototype.ipynb | 1777 ------------------------------------ rocof_test_data.py | 180 ++++ 3 files changed, 1983 insertions(+), 1777 deletions(-) create mode 100644 Phase Measurement Prototype.ipynb delete mode 100644 Phase measurement Prototype.ipynb create mode 100644 rocof_test_data.py diff --git a/Phase Measurement Prototype.ipynb b/Phase Measurement Prototype.ipynb new file mode 100644 index 0000000..df68f94 --- /dev/null +++ b/Phase Measurement Prototype.ipynb @@ -0,0 +1,1803 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 224, + "metadata": {}, + "outputs": [], + "source": [ + "import math\n", + "\n", + "import numpy as np\n", + "from scipy import signal, optimize\n", + "from matplotlib import pyplot as plt\n", + "\n", + "import rocof_test_data" + ] + }, + { + "cell_type": "code", + "execution_count": 225, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib notebook" + ] + }, + { + "cell_type": "code", + "execution_count": 255, + "metadata": { + "scrolled": false + }, + "outputs": [], + "source": [ + "fs = 1e3 # Hz\n", + "ff = 50 # Hz\n", + "duration = 60 # seconds\n", + "# test_data = rocof_test_data.sample_waveform(rocof_test_data.test_close_interharmonics_and_flicker(),\n", + "# duration=20,\n", + "# sampling_rate=fs,\n", + "# frequency=ff)[0]\n", + "# test_data = rocof_test_data.sample_waveform(rocof_test_data.gen_noise(fmin=10, amplitude=1),\n", + "# duration=20,\n", + "# sampling_rate=fs,\n", + "# frequency=ff)[0]\n", + "\n", + "\n", + "#gen = rocof_test_data.gen_noise(fmin=10, amplitude=1)\n", + "# gen = rocof_test_data.gen_noise(fmin=60, amplitude=0.2)\n", + "# gen = rocof_test_data.test_harmonics()\n", + "# gen = rocof_test_data.gen_interharmonic(*rocof_test_data.test_interharmonics)\n", + "# gen = rocof_test_data.test_amplitude_steps()\n", + "# gen = rocof_test_data.test_amplitude_and_phase_steps()\n", + "test_data = []\n", + "test_labels = [ fun.__name__.replace('test_', '') for fun in rocof_test_data.all_tests ]\n", + "for gen in rocof_test_data.all_tests:\n", + " test_data.append(rocof_test_data.sample_waveform(gen(),\n", + " duration=duration,\n", + " sampling_rate=fs,\n", + " frequency=ff)[0])\n", + "# d = 10 # seconds\n", + "# test_data = np.sin(2*np.pi * ff * np.linspace(0, d, int(d*fs)))" + ] + }, + { + "cell_type": "code", + "execution_count": 256, + "metadata": {}, + "outputs": [], + "source": [ + "analysis_periods = 10\n", + "window_len = fs * analysis_periods/ff\n", + "nfft_factor = 1\n", + "sigma = window_len/8 # samples\n", + "\n", + "ffts = []\n", + "for item in test_data:\n", + " f, t, Zxx = signal.stft(item,\n", + " fs = fs,\n", + " window=('gaussian', sigma),\n", + " nperseg = window_len,\n", + " nfft = window_len * nfft_factor)\n", + " #boundary = 'zeros')\n", + " ffts.append((f, t, Zxx))" + ] + }, + { + "cell_type": "code", + "execution_count": 257, + "metadata": { + "scrolled": false + }, + "outputs": [ + { + "data": { + "application/javascript": [ + "/* Put everything inside the global mpl namespace */\n", + "window.mpl = {};\n", + "\n", + "\n", + "mpl.get_websocket_type = function() {\n", + " if (typeof(WebSocket) !== 'undefined') {\n", + " return WebSocket;\n", + " } else if (typeof(MozWebSocket) !== 'undefined') {\n", + " return MozWebSocket;\n", + " } else {\n", + " alert('Your browser does not have WebSocket support.' +\n", + " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", + " 'Firefox 4 and 5 are also supported but you ' +\n", + " 'have to enable WebSockets in about:config.');\n", + " };\n", + "}\n", + "\n", + "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n", + " this.id = figure_id;\n", + "\n", + " this.ws = websocket;\n", + "\n", + " this.supports_binary = (this.ws.binaryType != undefined);\n", + "\n", + " if (!this.supports_binary) {\n", + " var warnings = document.getElementById(\"mpl-warnings\");\n", + " if (warnings) {\n", + " warnings.style.display = 'block';\n", + " warnings.textContent = (\n", + " \"This browser does not support binary websocket messages. \" +\n", + " \"Performance may be slow.\");\n", + " }\n", + " }\n", + "\n", + " this.imageObj = new Image();\n", + "\n", + " this.context = undefined;\n", + " this.message = undefined;\n", + " this.canvas = undefined;\n", + " this.rubberband_canvas = undefined;\n", + " this.rubberband_context = undefined;\n", + " this.format_dropdown = undefined;\n", + "\n", + " this.image_mode = 'full';\n", + "\n", + " this.root = $('
');\n", + " this._root_extra_style(this.root)\n", + " this.root.attr('style', 'display: inline-block');\n", + "\n", + " $(parent_element).append(this.root);\n", + "\n", + " this._init_header(this);\n", + " this._init_canvas(this);\n", + " this._init_toolbar(this);\n", + "\n", + " var fig = this;\n", + "\n", + " this.waiting = false;\n", + "\n", + " this.ws.onopen = function () {\n", + " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", + " fig.send_message(\"send_image_mode\", {});\n", + " if (mpl.ratio != 1) {\n", + " fig.send_message(\"set_dpi_ratio\", {'dpi_ratio': mpl.ratio});\n", + " }\n", + " fig.send_message(\"refresh\", {});\n", + " }\n", + "\n", + " this.imageObj.onload = function() {\n", + " if (fig.image_mode == 'full') {\n", + " // Full images could contain transparency (where diff images\n", + " // almost always do), so we need to clear the canvas so that\n", + " // there is no ghosting.\n", + " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", + " }\n", + " fig.context.drawImage(fig.imageObj, 0, 0);\n", + " };\n", + "\n", + " this.imageObj.onunload = function() {\n", + " fig.ws.close();\n", + " }\n", + "\n", + " this.ws.onmessage = this._make_on_message_function(this);\n", + "\n", + " this.ondownload = ondownload;\n", + "}\n", + "\n", + "mpl.figure.prototype._init_header = function() {\n", + " var titlebar = $(\n", + " '
');\n", + " var titletext = $(\n", + " '
');\n", + " titlebar.append(titletext)\n", + " this.root.append(titlebar);\n", + " this.header = titletext[0];\n", + "}\n", + "\n", + "\n", + "\n", + "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", + "\n", + "}\n", + "\n", + "\n", + "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", + "\n", + "}\n", + "\n", + "mpl.figure.prototype._init_canvas = function() {\n", + " var fig = this;\n", + "\n", + " var canvas_div = $('
');\n", + "\n", + " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", + "\n", + " function canvas_keyboard_event(event) {\n", + " return fig.key_event(event, event['data']);\n", + " }\n", + "\n", + " canvas_div.keydown('key_press', canvas_keyboard_event);\n", + " canvas_div.keyup('key_release', canvas_keyboard_event);\n", + " this.canvas_div = canvas_div\n", + " this._canvas_extra_style(canvas_div)\n", + " this.root.append(canvas_div);\n", + "\n", + " var canvas = $('');\n", + " canvas.addClass('mpl-canvas');\n", + " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", + "\n", + " this.canvas = canvas[0];\n", + " this.context = canvas[0].getContext(\"2d\");\n", + "\n", + " var backingStore = this.context.backingStorePixelRatio ||\n", + "\tthis.context.webkitBackingStorePixelRatio ||\n", + "\tthis.context.mozBackingStorePixelRatio ||\n", + "\tthis.context.msBackingStorePixelRatio ||\n", + "\tthis.context.oBackingStorePixelRatio ||\n", + "\tthis.context.backingStorePixelRatio || 1;\n", + "\n", + " mpl.ratio = (window.devicePixelRatio || 1) / backingStore;\n", + "\n", + " var rubberband = $('');\n", + " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", + "\n", + " var pass_mouse_events = true;\n", + "\n", + " canvas_div.resizable({\n", + " start: function(event, ui) {\n", + " pass_mouse_events = false;\n", + " },\n", + " resize: function(event, ui) {\n", + " fig.request_resize(ui.size.width, ui.size.height);\n", + " },\n", + " stop: function(event, ui) {\n", + " pass_mouse_events = true;\n", + " fig.request_resize(ui.size.width, ui.size.height);\n", + " },\n", + " });\n", + "\n", + " function mouse_event_fn(event) {\n", + " if (pass_mouse_events)\n", + " return fig.mouse_event(event, event['data']);\n", + " }\n", + "\n", + " rubberband.mousedown('button_press', mouse_event_fn);\n", + " rubberband.mouseup('button_release', mouse_event_fn);\n", + " // Throttle sequential mouse events to 1 every 20ms.\n", + " rubberband.mousemove('motion_notify', mouse_event_fn);\n", + "\n", + " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", + " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", + "\n", + " canvas_div.on(\"wheel\", function (event) {\n", + " event = event.originalEvent;\n", + " event['data'] = 'scroll'\n", + " if (event.deltaY < 0) {\n", + " event.step = 1;\n", + " } else {\n", + " event.step = -1;\n", + " }\n", + " mouse_event_fn(event);\n", + " });\n", + "\n", + " canvas_div.append(canvas);\n", + " canvas_div.append(rubberband);\n", + "\n", + " this.rubberband = rubberband;\n", + " this.rubberband_canvas = rubberband[0];\n", + " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", + " this.rubberband_context.strokeStyle = \"#000000\";\n", + "\n", + " this._resize_canvas = function(width, height) {\n", + " // Keep the size of the canvas, canvas container, and rubber band\n", + " // canvas in synch.\n", + " canvas_div.css('width', width)\n", + " canvas_div.css('height', height)\n", + "\n", + " canvas.attr('width', width * mpl.ratio);\n", + " canvas.attr('height', height * mpl.ratio);\n", + " canvas.attr('style', 'width: ' + width + 'px; height: ' + height + 'px;');\n", + "\n", + " rubberband.attr('width', width);\n", + " rubberband.attr('height', height);\n", + " }\n", + "\n", + " // Set the figure to an initial 600x600px, this will subsequently be updated\n", + " // upon first draw.\n", + " this._resize_canvas(600, 600);\n", + "\n", + " // Disable right mouse context menu.\n", + " $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", + " return false;\n", + " });\n", + "\n", + " function set_focus () {\n", + " canvas.focus();\n", + " canvas_div.focus();\n", + " }\n", + "\n", + " window.setTimeout(set_focus, 100);\n", + "}\n", + "\n", + "mpl.figure.prototype._init_toolbar = function() {\n", + " var fig = this;\n", + "\n", + " var nav_element = $('
')\n", + " nav_element.attr('style', 'width: 100%');\n", + " this.root.append(nav_element);\n", + "\n", + " // Define a callback function for later on.\n", + " function toolbar_event(event) {\n", + " return fig.toolbar_button_onclick(event['data']);\n", + " }\n", + " function toolbar_mouse_event(event) {\n", + " return fig.toolbar_button_onmouseover(event['data']);\n", + " }\n", + "\n", + " for(var toolbar_ind in mpl.toolbar_items) {\n", + " var name = mpl.toolbar_items[toolbar_ind][0];\n", + " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", + " var image = mpl.toolbar_items[toolbar_ind][2];\n", + " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", + "\n", + " if (!name) {\n", + " // put a spacer in here.\n", + " continue;\n", + " }\n", + " var button = $('');\n", + " button.click(method_name, toolbar_event);\n", + " button.mouseover(tooltip, toolbar_mouse_event);\n", + " nav_element.append(button);\n", + " }\n", + "\n", + " // Add the status bar.\n", + " var status_bar = $('');\n", + " nav_element.append(status_bar);\n", + " this.message = status_bar[0];\n", + "\n", + " // Add the close button to the window.\n", + " var buttongrp = $('
');\n", + " var button = $('');\n", + " button.click(function (evt) { fig.handle_close(fig, {}); } );\n", + " button.mouseover('Stop Interaction', toolbar_mouse_event);\n", + " buttongrp.append(button);\n", + " var titlebar = this.root.find($('.ui-dialog-titlebar'));\n", + " titlebar.prepend(buttongrp);\n", + "}\n", + "\n", + "mpl.figure.prototype._root_extra_style = function(el){\n", + " var fig = this\n", + " el.on(\"remove\", function(){\n", + "\tfig.close_ws(fig, {});\n", + " });\n", + "}\n", + "\n", + "mpl.figure.prototype._canvas_extra_style = function(el){\n", + " // this is important to make the div 'focusable\n", + " el.attr('tabindex', 0)\n", + " // reach out to IPython and tell the keyboard manager to turn it's self\n", + " // off when our div gets focus\n", + "\n", + " // location in version 3\n", + " if (IPython.notebook.keyboard_manager) {\n", + " IPython.notebook.keyboard_manager.register_events(el);\n", + " }\n", + " else {\n", + " // location in version 2\n", + " IPython.keyboard_manager.register_events(el);\n", + " }\n", + "\n", + "}\n", + "\n", + "mpl.figure.prototype._key_event_extra = function(event, name) {\n", + " var manager = IPython.notebook.keyboard_manager;\n", + " if (!manager)\n", + " manager = IPython.keyboard_manager;\n", + "\n", + " // Check for shift+enter\n", + " if (event.shiftKey && event.which == 13) {\n", + " this.canvas_div.blur();\n", + " event.shiftKey = false;\n", + " // Send a \"J\" for go to next cell\n", + " event.which = 74;\n", + " event.keyCode = 74;\n", + " manager.command_mode();\n", + " manager.handle_keydown(event);\n", + " }\n", + "}\n", + "\n", + "mpl.figure.prototype.handle_save = function(fig, msg) {\n", + " fig.ondownload(fig, null);\n", + "}\n", + "\n", + "\n", + "mpl.find_output_cell = function(html_output) {\n", + " // Return the cell and output element which can be found *uniquely* in the notebook.\n", + " // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n", + " // IPython event is triggered only after the cells have been serialised, which for\n", + " // our purposes (turning an active figure into a static one), is too late.\n", + " var cells = IPython.notebook.get_cells();\n", + " var ncells = cells.length;\n", + " for (var i=0; i= 3 moved mimebundle to data attribute of output\n", + " data = data.data;\n", + " }\n", + " if (data['text/html'] == html_output) {\n", + " return [cell, data, j];\n", + " }\n", + " }\n", + " }\n", + " }\n", + "}\n", + "\n", + "// Register the function which deals with the matplotlib target/channel.\n", + "// The kernel may be null if the page has been refreshed.\n", + "if (IPython.notebook.kernel != null) {\n", + " IPython.notebook.kernel.comm_manager.register_target('matplotlib', mpl.mpl_figure_comm);\n", + "}\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "fig, axs = plt.subplots(len(test_data), figsize=(8, 20), sharex=True)\n", + "fig.tight_layout(pad=2.2, h_pad=0, w_pad=1)\n", + "\n", + "for fft, ax, label in zip(ffts, axs.flatten(), test_labels):\n", + " f, f_t, Zxx = fft\n", + " \n", + " n_f, n_t = Zxx.shape\n", + " # start, stop = 180, 220\n", + " # start, stop = 90, 110\n", + " # start, stop = 15, 35\n", + " # bounds_f = slice(start // 4 * nfft_factor, stop // 4 * nfft_factor)\n", + " f_min, f_max = 30, 70 # Hz\n", + " bounds_f = slice(np.argmax(f > f_min), np.argmin(f < f_max))\n", + " \n", + "\n", + " f_mean = np.zeros(Zxx.shape[1])\n", + " for t in range(1, Zxx.shape[1] - 1):\n", + " frame_f = f[bounds_f]\n", + " frame_step = frame_f[1] - frame_f[0]\n", + " time_step = f_t[1] - f_t[0]\n", + " #if t == 10:\n", + " # axs[-1].plot(frame_f, frame_Z)\n", + " frame_Z = np.abs(Zxx[bounds_f, t])\n", + " # frame_f = f[180:220]\n", + " # frame_Z = np.abs(Zxx[180:220, 40])\n", + " # frame_f = f[15:35]\n", + " # frame_Z = np.abs(Zxx[15:35, 40])\n", + " # plt.plot(frame_f, frame_Z)\n", + "\n", + " # peak_f = frame_f[np.argmax(frame)]\n", + " # plt.axvline(peak_f, color='red')\n", + "\n", + "# def gauss(x, *p):\n", + "# A, mu, sigma, o = p\n", + "# return A*np.exp(-(x-mu)**2/(2.*sigma**2)) + o\n", + " \n", + " def gauss(x, *p):\n", + " A, mu, sigma = p\n", + " return A*np.exp(-(x-mu)**2/(2.*sigma**2))\n", + " \n", + " f_start = frame_f[np.argmax(frame_Z)]\n", + " A_start = np.max(frame_Z)\n", + " p0 = [A_start, f_start, 1.]\n", + " try:\n", + " coeff, var = optimize.curve_fit(gauss, frame_f, frame_Z, p0=p0)\n", + " # plt.plot(frame_f, gauss(frame_f, *coeff))\n", + " #print(coeff)\n", + " A, mu, sigma, *_ = coeff\n", + " f_mean[t] = mu\n", + " except RuntimeError:\n", + " f_mean[t] = np.nan\n", + " ax.plot(f_t[1:-1], f_mean[1:-1])\n", + " \n", + "# b, a = signal.butter(3,\n", + "# 1/5, # Hz\n", + "# btype='lowpass',\n", + "# fs=1/time_step)\n", + "# filtered = signal.lfilter(b, a, f_mean[1:-1], axis=0)\n", + "# ax.plot(f_t[1:-1], filtered)\n", + " \n", + " ax.set_title(label, pad=-20)\n", + " ax.set_ylabel('f [Hz]')\n", + " ax.grid()\n", + " if not label in ['off_frequency', 'sweep_phase_steps']:\n", + " ax.set_ylim([49.90, 50.10])\n", + " var = np.var(f_mean[1:-1])\n", + " ax.text(0.5, 0.1, f'σ²={var * 1e3:.3g} mHz²', transform=ax.transAxes, ha='center')\n", + " ax.text(0.5, 0.25, f'σ={np.sqrt(var) * 1e3:.3g} mHz', transform=ax.transAxes, ha='center')\n", + "# ax.text(0.5, 0.2, f'filt. σ²={np.var(filtered) * 1e3:.3g} mHz', transform=ax.transAxes, ha='center')\n", + " else:\n", + " f_min, f_max = min(f_mean[1:-1]), max(f_mean[1:-1])\n", + " delta = f_max - f_min\n", + " ax.set_ylim(f_min - delta * 0.1, f_max + delta * 0.3)\n", + " \n", + "ax.set_xlabel('simulation time t [s]')\n", + "None" + ] + } + ], + "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.7.5" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/Phase measurement Prototype.ipynb b/Phase measurement Prototype.ipynb deleted file mode 100644 index 498f274..0000000 --- a/Phase measurement Prototype.ipynb +++ /dev/null @@ -1,1777 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": 7, - "metadata": {}, - "outputs": [], - "source": [ - "import math\n", - "\n", - "import numpy as np\n", - "from scipy import signal, optimize\n", - "from matplotlib import pyplot as plt\n", - "\n", - "import rocof_test_data" - ] - }, - { - "cell_type": "code", - "execution_count": 8, - "metadata": {}, - "outputs": [], - "source": [ - "%matplotlib notebook" - ] - }, - { - "cell_type": "code", - "execution_count": 9, - "metadata": { - "scrolled": false - }, - "outputs": [], - "source": [ - "fs = 1e3 # Hz\n", - "ff = 50 # Hz\n", - "duration = 60 # seconds\n", - "# test_data = rocof_test_data.sample_waveform(rocof_test_data.test_close_interharmonics_and_flicker(),\n", - "# duration=20,\n", - "# sampling_rate=fs,\n", - "# frequency=ff)[0]\n", - "# test_data = rocof_test_data.sample_waveform(rocof_test_data.gen_noise(fmin=10, amplitude=1),\n", - "# duration=20,\n", - "# sampling_rate=fs,\n", - "# frequency=ff)[0]\n", - "\n", - "\n", - "#gen = rocof_test_data.gen_noise(fmin=10, amplitude=1)\n", - "# gen = rocof_test_data.gen_noise(fmin=60, amplitude=0.2)\n", - "# gen = rocof_test_data.test_harmonics()\n", - "# gen = rocof_test_data.gen_interharmonic(*rocof_test_data.test_interharmonics)\n", - "# gen = rocof_test_data.test_amplitude_steps()\n", - "# gen = rocof_test_data.test_amplitude_and_phase_steps()\n", - "test_data = []\n", - "test_labels = [ fun.__name__.replace('test_', '') for fun in rocof_test_data.all_tests ]\n", - "for gen in rocof_test_data.all_tests:\n", - " test_data.append(rocof_test_data.sample_waveform(gen(),\n", - " duration=duration,\n", - " sampling_rate=fs,\n", - " frequency=ff)[0])\n", - "# d = 10 # seconds\n", - "# test_data = np.sin(2*np.pi * ff * np.linspace(0, d, int(d*fs)))" - ] - }, - { - "cell_type": "code", - "execution_count": 10, - "metadata": {}, - "outputs": [], - "source": [ - "analysis_periods = 50\n", - "window_len = fs * analysis_periods/ff\n", - "nfft_factor = 1\n", - "\n", - "ffts = []\n", - "for item in test_data:\n", - " f, t, Zxx = signal.stft(item,\n", - " fs = fs,\n", - " window='blackman',\n", - " nperseg = window_len,\n", - " nfft = window_len * nfft_factor)\n", - " #boundary = 'zeros')\n", - " ffts.append((f, t, Zxx))" - ] - }, - { - "cell_type": "code", - "execution_count": 11, - "metadata": { - "scrolled": false - }, - "outputs": [ - { - "data": { - "application/javascript": [ - "/* Put everything inside the global mpl namespace */\n", - "window.mpl = {};\n", - "\n", - "\n", - "mpl.get_websocket_type = function() {\n", - " if (typeof(WebSocket) !== 'undefined') {\n", - " return WebSocket;\n", - " } else if (typeof(MozWebSocket) !== 'undefined') {\n", - " return MozWebSocket;\n", - " } else {\n", - " alert('Your browser does not have WebSocket support.' +\n", - " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", - " 'Firefox 4 and 5 are also supported but you ' +\n", - " 'have to enable WebSockets in about:config.');\n", - " };\n", - "}\n", - "\n", - "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n", - " this.id = figure_id;\n", - "\n", - " this.ws = websocket;\n", - "\n", - " this.supports_binary = (this.ws.binaryType != undefined);\n", - "\n", - " if (!this.supports_binary) {\n", - " var warnings = document.getElementById(\"mpl-warnings\");\n", - " if (warnings) {\n", - " warnings.style.display = 'block';\n", - " warnings.textContent = (\n", - " \"This browser does not support binary websocket messages. \" +\n", - " \"Performance may be slow.\");\n", - " }\n", - " }\n", - "\n", - " this.imageObj = new Image();\n", - "\n", - " this.context = undefined;\n", - " this.message = undefined;\n", - " this.canvas = undefined;\n", - " this.rubberband_canvas = undefined;\n", - " this.rubberband_context = undefined;\n", - " this.format_dropdown = undefined;\n", - "\n", - " this.image_mode = 'full';\n", - "\n", - " this.root = $('
');\n", - " this._root_extra_style(this.root)\n", - " this.root.attr('style', 'display: inline-block');\n", - "\n", - " $(parent_element).append(this.root);\n", - "\n", - " this._init_header(this);\n", - " this._init_canvas(this);\n", - " this._init_toolbar(this);\n", - "\n", - " var fig = this;\n", - "\n", - " this.waiting = false;\n", - "\n", - " this.ws.onopen = function () {\n", - " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", - " fig.send_message(\"send_image_mode\", {});\n", - " if (mpl.ratio != 1) {\n", - " fig.send_message(\"set_dpi_ratio\", {'dpi_ratio': mpl.ratio});\n", - " }\n", - " fig.send_message(\"refresh\", {});\n", - " }\n", - "\n", - " this.imageObj.onload = function() {\n", - " if (fig.image_mode == 'full') {\n", - " // Full images could contain transparency (where diff images\n", - " // almost always do), so we need to clear the canvas so that\n", - " // there is no ghosting.\n", - " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", - " }\n", - " fig.context.drawImage(fig.imageObj, 0, 0);\n", - " };\n", - "\n", - " this.imageObj.onunload = function() {\n", - " fig.ws.close();\n", - " }\n", - "\n", - " this.ws.onmessage = this._make_on_message_function(this);\n", - "\n", - " this.ondownload = ondownload;\n", - "}\n", - "\n", - "mpl.figure.prototype._init_header = function() {\n", - " var titlebar = $(\n", - " '
');\n", - " var titletext = $(\n", - " '
');\n", - " titlebar.append(titletext)\n", - " this.root.append(titlebar);\n", - " this.header = titletext[0];\n", - "}\n", - "\n", - "\n", - "\n", - "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", - "\n", - "}\n", - "\n", - "\n", - "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", - "\n", - "}\n", - "\n", - "mpl.figure.prototype._init_canvas = function() {\n", - " var fig = this;\n", - "\n", - " var canvas_div = $('
');\n", - "\n", - " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", - "\n", - " function canvas_keyboard_event(event) {\n", - " return fig.key_event(event, event['data']);\n", - " }\n", - "\n", - " canvas_div.keydown('key_press', canvas_keyboard_event);\n", - " canvas_div.keyup('key_release', canvas_keyboard_event);\n", - " this.canvas_div = canvas_div\n", - " this._canvas_extra_style(canvas_div)\n", - " this.root.append(canvas_div);\n", - "\n", - " var canvas = $('');\n", - " canvas.addClass('mpl-canvas');\n", - " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", - "\n", - " this.canvas = canvas[0];\n", - " this.context = canvas[0].getContext(\"2d\");\n", - "\n", - " var backingStore = this.context.backingStorePixelRatio ||\n", - "\tthis.context.webkitBackingStorePixelRatio ||\n", - "\tthis.context.mozBackingStorePixelRatio ||\n", - "\tthis.context.msBackingStorePixelRatio ||\n", - "\tthis.context.oBackingStorePixelRatio ||\n", - "\tthis.context.backingStorePixelRatio || 1;\n", - "\n", - " mpl.ratio = (window.devicePixelRatio || 1) / backingStore;\n", - "\n", - " var rubberband = $('');\n", - " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", - "\n", - " var pass_mouse_events = true;\n", - "\n", - " canvas_div.resizable({\n", - " start: function(event, ui) {\n", - " pass_mouse_events = false;\n", - " },\n", - " resize: function(event, ui) {\n", - " fig.request_resize(ui.size.width, ui.size.height);\n", - " },\n", - " stop: function(event, ui) {\n", - " pass_mouse_events = true;\n", - " fig.request_resize(ui.size.width, ui.size.height);\n", - " },\n", - " });\n", - "\n", - " function mouse_event_fn(event) {\n", - " if (pass_mouse_events)\n", - " return fig.mouse_event(event, event['data']);\n", - " }\n", - "\n", - " rubberband.mousedown('button_press', mouse_event_fn);\n", - " rubberband.mouseup('button_release', mouse_event_fn);\n", - " // Throttle sequential mouse events to 1 every 20ms.\n", - " rubberband.mousemove('motion_notify', mouse_event_fn);\n", - "\n", - " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", - " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", - "\n", - " canvas_div.on(\"wheel\", function (event) {\n", - " event = event.originalEvent;\n", - " event['data'] = 'scroll'\n", - " if (event.deltaY < 0) {\n", - " event.step = 1;\n", - " } else {\n", - " event.step = -1;\n", - " }\n", - " mouse_event_fn(event);\n", - " });\n", - "\n", - " canvas_div.append(canvas);\n", - " canvas_div.append(rubberband);\n", - "\n", - " this.rubberband = rubberband;\n", - " this.rubberband_canvas = rubberband[0];\n", - " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", - " this.rubberband_context.strokeStyle = \"#000000\";\n", - "\n", - " this._resize_canvas = function(width, height) {\n", - " // Keep the size of the canvas, canvas container, and rubber band\n", - " // canvas in synch.\n", - " canvas_div.css('width', width)\n", - " canvas_div.css('height', height)\n", - "\n", - " canvas.attr('width', width * mpl.ratio);\n", - " canvas.attr('height', height * mpl.ratio);\n", - " canvas.attr('style', 'width: ' + width + 'px; height: ' + height + 'px;');\n", - "\n", - " rubberband.attr('width', width);\n", - " rubberband.attr('height', height);\n", - " }\n", - "\n", - " // Set the figure to an initial 600x600px, this will subsequently be updated\n", - " // upon first draw.\n", - " this._resize_canvas(600, 600);\n", - "\n", - " // Disable right mouse context menu.\n", - " $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", - " return false;\n", - " });\n", - "\n", - " function set_focus () {\n", - " canvas.focus();\n", - " canvas_div.focus();\n", - " }\n", - "\n", - " window.setTimeout(set_focus, 100);\n", - "}\n", - "\n", - "mpl.figure.prototype._init_toolbar = function() {\n", - " var fig = this;\n", - "\n", - " var nav_element = $('
')\n", - " nav_element.attr('style', 'width: 100%');\n", - " this.root.append(nav_element);\n", - "\n", - " // Define a callback function for later on.\n", - " function toolbar_event(event) {\n", - " return fig.toolbar_button_onclick(event['data']);\n", - " }\n", - " function toolbar_mouse_event(event) {\n", - " return fig.toolbar_button_onmouseover(event['data']);\n", - " }\n", - "\n", - " for(var toolbar_ind in mpl.toolbar_items) {\n", - " var name = mpl.toolbar_items[toolbar_ind][0];\n", - " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", - " var image = mpl.toolbar_items[toolbar_ind][2];\n", - " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", - "\n", - " if (!name) {\n", - " // put a spacer in here.\n", - " continue;\n", - " }\n", - " var button = $('');\n", - " button.click(method_name, toolbar_event);\n", - " button.mouseover(tooltip, toolbar_mouse_event);\n", - " nav_element.append(button);\n", - " }\n", - "\n", - " // Add the status bar.\n", - " var status_bar = $('');\n", - " nav_element.append(status_bar);\n", - " this.message = status_bar[0];\n", - "\n", - " // Add the close button to the window.\n", - " var buttongrp = $('
');\n", - " var button = $('');\n", - " button.click(function (evt) { fig.handle_close(fig, {}); } );\n", - " button.mouseover('Stop Interaction', toolbar_mouse_event);\n", - " buttongrp.append(button);\n", - " var titlebar = this.root.find($('.ui-dialog-titlebar'));\n", - " titlebar.prepend(buttongrp);\n", - "}\n", - "\n", - "mpl.figure.prototype._root_extra_style = function(el){\n", - " var fig = this\n", - " el.on(\"remove\", function(){\n", - "\tfig.close_ws(fig, {});\n", - " });\n", - "}\n", - "\n", - "mpl.figure.prototype._canvas_extra_style = function(el){\n", - " // this is important to make the div 'focusable\n", - " el.attr('tabindex', 0)\n", - " // reach out to IPython and tell the keyboard manager to turn it's self\n", - " // off when our div gets focus\n", - "\n", - " // location in version 3\n", - " if (IPython.notebook.keyboard_manager) {\n", - " IPython.notebook.keyboard_manager.register_events(el);\n", - " }\n", - " else {\n", - " // location in version 2\n", - " IPython.keyboard_manager.register_events(el);\n", - " }\n", - "\n", - "}\n", - "\n", - "mpl.figure.prototype._key_event_extra = function(event, name) {\n", - " var manager = IPython.notebook.keyboard_manager;\n", - " if (!manager)\n", - " manager = IPython.keyboard_manager;\n", - "\n", - " // Check for shift+enter\n", - " if (event.shiftKey && event.which == 13) {\n", - " this.canvas_div.blur();\n", - " event.shiftKey = false;\n", - " // Send a \"J\" for go to next cell\n", - " event.which = 74;\n", - " event.keyCode = 74;\n", - " manager.command_mode();\n", - " manager.handle_keydown(event);\n", - " }\n", - "}\n", - "\n", - "mpl.figure.prototype.handle_save = function(fig, msg) {\n", - " fig.ondownload(fig, null);\n", - "}\n", - "\n", - "\n", - "mpl.find_output_cell = function(html_output) {\n", - " // Return the cell and output element which can be found *uniquely* in the notebook.\n", - " // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n", - " // IPython event is triggered only after the cells have been serialised, which for\n", - " // our purposes (turning an active figure into a static one), is too late.\n", - " var cells = IPython.notebook.get_cells();\n", - " var ncells = cells.length;\n", - " for (var i=0; i= 3 moved mimebundle to data attribute of output\n", - " data = data.data;\n", - " }\n", - " if (data['text/html'] == html_output) {\n", - " return [cell, data, j];\n", - " }\n", - " }\n", - " }\n", - " }\n", - "}\n", - "\n", - "// Register the function which deals with the matplotlib target/channel.\n", - "// The kernel may be null if the page has been refreshed.\n", - "if (IPython.notebook.kernel != null) {\n", - " IPython.notebook.kernel.comm_manager.register_target('matplotlib', mpl.mpl_figure_comm);\n", - "}\n" - ], - "text/plain": [ - "" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/html": [ - "" - ], - "text/plain": [ - "" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], - "source": [ - "fig, axs = plt.subplots(len(test_data), figsize=(8, 20), sharex=True)\n", - "fig.tight_layout(pad=2.2, h_pad=0, w_pad=1)\n", - "\n", - "for fft, ax, label in zip(ffts, axs.flatten(), test_labels):\n", - " f, f_t, Zxx = fft\n", - " \n", - " n_f, n_t = Zxx.shape\n", - " start, stop = 180, 220\n", - " # start, stop = 90, 110\n", - " # start, stop = 15, 35\n", - " bounds_f = slice(start // 4 * nfft_factor, stop // 4 * nfft_factor)\n", - " \n", - "\n", - " f_mean = np.zeros(Zxx.shape[1])\n", - " for t in range(1, Zxx.shape[1] - 1):\n", - " frame_f = f[bounds_f]\n", - " frame_Z = np.abs(Zxx[bounds_f, t])\n", - " # frame_f = f[180:220]\n", - " # frame_Z = np.abs(Zxx[180:220, 40])\n", - " # frame_f = f[15:35]\n", - " # frame_Z = np.abs(Zxx[15:35, 40])\n", - " # plt.plot(frame_f, frame_Z)\n", - "\n", - " # peak_f = frame_f[np.argmax(frame)]\n", - " # plt.axvline(peak_f, color='red')\n", - "\n", - " def gauss(x, *p):\n", - " A, mu, sigma, o = p\n", - " return A*np.exp(-(x-mu)**2/(2.*sigma**2)) + o\n", - " \n", - " f_start = frame_f[np.argmax(frame_Z)]\n", - " p0 = [.5, f_start, 1., 0]\n", - " coeff, var = optimize.curve_fit(gauss, frame_f, frame_Z, p0=p0)\n", - " # plt.plot(frame_f, gauss(frame_f, *coeff))\n", - " #print(coeff)\n", - " A, mu, sigma, o = coeff\n", - " f_mean[t] = mu\n", - " ax.plot(f_t[1:-1], f_mean[1:-1])\n", - " ax.set_title(label, pad=-20)\n", - " ax.set_ylabel('f [Hz]')\n", - " ax.grid()\n", - " if not label in ['off_frequency', 'sweep_phase_steps']:\n", - " ax.set_ylim([49.90, 50.10])\n", - " ax.text(0.5, 0.1, f'σ²={np.var(f_mean[1:-1]) * 1e3:.3g} mHz', transform=ax.transAxes, ha='center')\n", - " else:\n", - " f_min, f_max = min(f_mean[1:-1]), max(f_mean[1:-1])\n", - " delta = f_max - f_min\n", - " ax.set_ylim(f_min - delta * 0.1, f_max + delta * 0.3)\n", - " \n", - "ax.set_xlabel('simulation time t [s]')\n", - "None" - ] - } - ], - "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.7.5" - } - }, - "nbformat": 4, - "nbformat_minor": 2 -} diff --git a/rocof_test_data.py b/rocof_test_data.py new file mode 100644 index 0000000..91bee95 --- /dev/null +++ b/rocof_test_data.py @@ -0,0 +1,180 @@ +#!/usr/bin/env python +# coding: utf-8 + +# # ROCOF test waveform library +# +# This is a re-implementation of the ROCOF test waveforms described in https://zenodo.org/record/3559798 +# +# **This file is exported as a python module and loaded from other notebooks here, so please make sure to re-export when changing it.** + +# In[ ]: + + +import math +import itertools + +import numpy as np +from scipy import signal +from matplotlib import pyplot as plt + + +# In[ ]: + + +get_ipython().run_line_magic('matplotlib', 'notebook') + + +# In[ ]: + + +def sample_waveform(generator, duration:"s"=10, sampling_rate:"sp/s"=10000, frequency:"Hz"=50): + samples = int(duration*sampling_rate) + phases = np.linspace(0, 2*np.pi, 6, endpoint=False) + omega_t = np.linspace(phases, phases + 2*np.pi*duration*frequency, samples) + fundamental = np.sin(omega_t) + return generator(omega_t, fundamental, sampling_rate=sampling_rate, duration=duration, frequency=frequency).swapaxes(0, 1) + + +# In[ ]: + + +def gen_harmonics(amplitudes, phases=[]): + return lambda omega_t, fundamental, **_: fundamental + np.sum([ + a*np.sin((p if p else 0) + i*omega_t) + for i, (a, p) in enumerate(itertools.zip_longest(amplitudes, phases), start=2) + ], axis=0) + +def test_harmonics(): + return gen_harmonics([0.02, 0.05, 0.01, 0.06, 0.005, 0.05, 0.005, 0.015, 0.005, 0.035, 0.005, 0.003]) + + +# In[ ]: + + +def gen_interharmonic(amplitudes, ih=[], ih_phase=[]): + def gen(omega_t, fundamental, **_): + return fundamental + np.sum([ + a*np.sin(omega_t * ih + (p if p else 0)) + for a, ih, p in itertools.zip_longest(amplitudes, ih, ih_phase) + ], axis=0) + return gen + +def test_interharmonics(): + return gen_interharmonic([0.1], [15.01401], [np.pi]) + + +# In[ ]: + + +def gen_noise(amplitude=0.2, fmax:'Hz'=4.9e3, fmin:'Hz'=100, filter_order=6): + def gen(omega_t, fundamental, sampling_rate, **_): + noise = np.random.normal(0, amplitude, fundamental.shape) + b, a = signal.butter(filter_order, + [fmin, min(fmax, sampling_rate//2-1)], + btype='bandpass', + fs=sampling_rate) + return fundamental + signal.lfilter(b, a, noise, axis=0) + return gen + +def test_noise(): + return gen_noise() + +def test_noise_loud(): + return gen_noise(amplitude=0.5, fmin=10) + + +# In[406]: + + +def gen_steps(size_amplitude=0.1, size_phase=0.1*np.pi, steps_per_sec=1): + def gen(omega_t, fundamental, duration, **_): + n = int(steps_per_sec * duration) + indices = np.random.randint(0, len(omega_t), n) + amplitudes = np.random.normal(1, size_amplitude, (n, 6)) + phases = np.random.normal(0, size_phase, (n, 6)) + amplitude = np.ones(omega_t.shape) + for start, end, a, p in zip(indices, indices[1:], amplitudes, phases): + omega_t[start:end] += p + amplitude[start:end] = a + return amplitude*np.sin(omega_t) + return gen + +def test_amplitude_steps(): + return gen_steps(size_amplitude=0.4, size_phase=0) + +def test_phase_steps(): + return gen_steps(size_amplitude=0, size_phase=0.1) + +def test_amplitude_and_phase_steps(): + return gen_steps(size_amplitude=0.2, size_phase=0.07) + + +# In[418]: + + +def step_gen(shape, stdev, duration, steps_per_sec=1.0, mean=0.0): + samples, channels = shape + n = int(steps_per_sec * duration) + indices = np.random.randint(0, samples, n) + phases = np.random.normal(mean, stdev, (n, 6)) + amplitude = np.ones((samples, channels)) + out = np.zeros(shape) + for start, end, a in zip(indices, indices[1:], amplitude): + out[start:end] = a + return out + +def gen_chirp(fmin, fmax, period, dwell_time=1.0, amplitude=None, phase_steps=None): + def gen(omega_t, fundamental, sampling_rate, duration, **_): + samples = int(duration*sampling_rate) + phases = np.linspace(0, 2*np.pi, 6, endpoint=False) + + c = (fmax-fmin)/period + t = np.linspace(0, duration, samples) + + x = np.repeat(np.reshape(2*np.pi*fmin*t, (-1,1)), 6, axis=1) + data = (phases + x)[:int(sampling_rate*dwell_time)] + current_phase = 2*np.pi*fmin*dwell_time + direction = 'up' + + for idx in range(int(dwell_time*sampling_rate), samples, int(2*period*sampling_rate)): + t1 = np.linspace(0, period, int(period*sampling_rate)) + t2 = np.linspace(0, period, int(period*sampling_rate)) + chirp_phase = np.hstack(( + 2*np.pi*(c/2 * t1**2 + fmin * t1), + 2*np.pi*(-c/2 * t2**2 + fmax * t2 - (c/2 * period**2 + fmin * period)) + )) + chirp_phase = np.repeat(np.reshape(chirp_phase, (-1, 1)), 6, axis=1) + new = phases + chirp_phase + current_phase + current_phase = chirp_phase[-1] + data = np.vstack((data, new)) + + data = data[:len(fundamental)] + + if phase_steps: + (step_amplitude, steps_per_sec) = phase_steps + steps = step_gen(data.shape, step_amplitude, duration, steps_per_sec) + data += steps + + if amplitude is None: + return np.sin(data) + else: + return fundamental + amplitude*np.sin(data) + return gen + +def test_close_interharmonics_and_flicker(): + return gen_chirp(90.0, 150.0, 10, 1, amplitude=0.1) + +def test_off_frequency(): +# return gen_chirp(48.0, 52.0, 0.25, 1) + return gen_chirp(48.0, 52.0, 10, 1) + +def test_sweep_phase_steps(): + return gen_chirp(48.0, 52.0, 10, 1, phase_steps=(0.1, 1)) +# return gen_chirp(48.0, 52.0, 0.25, 1, phase_steps=(0.1, 1)) + + +# In[ ]: + + +all_tests = [test_harmonics, test_interharmonics, test_noise, test_noise_loud, test_amplitude_steps, test_phase_steps, test_amplitude_and_phase_steps, test_close_interharmonics_and_flicker, test_off_frequency, test_sweep_phase_steps] + -- cgit