{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Get some ASC data - Calculate Sensing Matrix \n", "### also make the radar plots" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Library Imports and Python parameter settings\n", "%matplotlib notebook\n", "#from __future__ import division\n", "import nds2\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import scipy.signal as sig\n", "\n", "import h5py\n", "import deepdish\n", "from astropy.time import Time\n", "import time\n", "import timeit\n", "\n", "# only works in python 2\n", "# import sfft\n", "\n", "plt.style.use('seaborn-whitegrid')\n", "plt.style.use('seaborn-paper')\n", "# Update the matplotlib configuration parameters:\n", "plt.rcParams.update({'text.usetex': False,\n", " 'lines.linewidth': 3,\n", " 'font.family': 'sans-serif',\n", " 'font.serif': 'Helvetica',\n", " 'font.size': 14,\n", " 'xtick.labelsize': 'x-large',\n", " 'ytick.labelsize': 'x-large',\n", " 'axes.labelsize': 'x-large',\n", " 'axes.titlesize': 'x-large',\n", " 'axes.grid': True,\n", " 'grid.alpha': 0.53,\n", " 'lines.markersize': 12,\n", " 'legend.borderpad': 0.2,\n", " 'legend.fancybox': True,\n", " 'legend.fontsize': 'small',\n", " 'legend.framealpha': 0.7,\n", " 'legend.handletextpad': 0.1,\n", " 'legend.labelspacing': 0.2,\n", " 'legend.loc': 'best',\n", " 'figure.figsize': (12,8),\n", " 'savefig.dpi': 100,\n", " 'pdf.compression': 9})\n", "plt.style.use('fivethirtyeight')\n", "\n", "debugme = 1\n", "save_data = 0\n", "save_plots = 0" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### How does this all work?\n", "* setup the NDS connection and \n", "* construct a list of channels and \n", "* get the data" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Setup connection to the NDS\n", "# needs Kerberos to make authentic conn to NDS2 servers at LLO, LHO, or CIT\n", "# *** takes ~5-10 minutes to download the whole channel list on the first connect of the day ***\n", "print('Making NDS3 connection...')\n", "#conn = nds2.connection('nds.ligo-la.caltech.edu', 31200) # LLO data from outside\n", "#conn = nds2.connection('nds.ligo-wa.caltech.edu', 31200) # LHO data from outside\n", "conn = nds2.connection('nds.ligo.caltech.edu', 31200) # LIGO data from outside\n", "#conn = nds2.connection('l1nds1', 8088) # LLO data from the inside\n", "\n", "# Setup start and stop times\n", "#times = '2018-12-20 10:09:00'\n", "times = '2018-12-23 9:07:00'\n", "t = Time(times, format='iso', scale='utc')\n", "t_start = int(t.gps)\n", "dur = 30 * 60\n", "\n", "# channel names\n", "chan_pref = 'H1:LSC-'\n", "chan_suff = '_ERR_DQ'\n", "ports = ['REFL','POP']\n", "freqs = ['9','45']\n", "dets = ['A']\n", "phases = ['I','Q']\n", "chans=[]\n", "for port in ports:\n", " for freq in freqs:\n", " for det in dets:\n", " for pha in phases:\n", " chans.append(chan_pref + port + '_' + det + '_' + 'RF' + freq + '_' + pha + chan_suff)\n", "\n", "channels = chans\n", "if debugme:\n", " for k in channels:\n", " print(k)\n", "\n", "# ~~~~ Get the data ~~~~~\n", "tic = timeit.default_timer()\n", "if debugme:\n", " print('Fetching data...')\n", "data = conn.fetch(t_start, t_start + dur, channels)\n", "toc = timeit.default_timer() - tic\n", "if debugme:\n", " print('Got data: Elapsed time = ' + str(round(toc,1)) + ' s')\n", " \n", "if save_data:\n", " print(\"Save data to a file.\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# some normalization & remove the filter transient\n", "t_cut = 1\n", "f_s = int(data[0].sample_rate)\n", "# Do some signal conditioning\n", "misos = sig.cheby1(8, 0.1, [100,200], btype='bandpass', fs=f_s, output='sos')\n", "vdata = []\n", "\n", "g_norm = [1e5, # PRC1\n", " 50, # PRC2\n", " 100, # SRC1\n", " 1e5, # SRC2\n", " 1, # CHARD\n", " 1e5] # CSOFT\n", "\n", "chan_idx = range(len(channels))\n", "for k in chan_idx:\n", " y = data[k].data# * g_norm[k]\n", " y = sig.sosfiltfilt(misos, y)\n", " y = y[t_cut*f_s:-t_cut*f_s]\n", " y = sig.detrend(y)\n", " vdata.append(y)\n", "\n", "y = np.asarray(vdata)\n", "\n", "# make ASD from time series data, apply inverse whitening to spectrum\n", "nFFT = f_s * 64\n", "\n", "#f, Phh = sig.welch(h1bs, fs=f_s, nperseg=nFFT, noverlap = nFFT/2, window=sig.get_window('hann', nFFT))\n", "f, Pyy = sig.welch(y, fs=f_s, nperseg=nFFT, noverlap = nFFT/2, \n", " window=sig.get_window('hann', nFFT), average='median')\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Make PSDs of the filtered data and see the peaks" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "fig = plt.figure(figsize=(20,10))\n", "\n", "Ayy = np.sqrt(Pyy)\n", "#plt.semilogy(f, Ahh, linestyle='-', lw=2, alpha=0.90, color='Green',\n", "# label=\"H1\")\n", "for k in range(len(Ayy)):\n", " plt.semilogy(f, Ayy[k,:], \n", " linestyle='-', lw=2, alpha=0.7, label = chans[chan_idx[k]])\n", "\n", "plt.xlabel(r'Frequency [Hz]', fontsize=18)\n", "plt.ylabel(r'Displacment [cts/$\\sqrt{Hz}$]', fontsize=18)\n", "\n", "#plt.title(\"IR Transmission during ALS Scan\");\n", "plt.grid(b=True, which=\"major\")\n", "plt.grid(b=True, which=\"minor\",alpha=0.21)\n", "plt.axis('tight')\n", "plt.xlim((121, 179))\n", "plt.ylim((1e-7, 31))\n", "\n", "plt.legend()\n", "plt.savefig('Data/' + 'MahLines.pdf', bbox_inches='tight')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Demodulate the time series data" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "```\n", "PRM 152 Hz\n", "SRM 138 Hz\n", "BS 128 Hz\n", "```" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "freqs = [152, 138, 128]\n", "phi = [45, 90, -10]\n", "exc = ['PRM', 'SRM', 'BS']\n", "\n", "a,b = y.shape\n", "\n", "down_fac = 256\n", "lpsos = sig.butter(2, 0.003, btype='lowpass', fs=int(f_s/down_fac), output='sos')\n", "\n", "t = np.linspace(0, b/f_s, b)\n", "IQ = np.zeros((a,len(freqs),int(b/down_fac)), dtype=complex) # {chans, freqs, t}\n", "\n", "for j in range(len(y)):\n", " for k in range(len(freqs)):\n", " ello = np.exp(1j * (2 * np.pi * freqs[k] * t + np.pi/180*phi[k])) # the local osc\n", " x = y[j] * ello # the mixer\n", " x = sig.decimate(x, down_fac, ftype='fir', zero_phase=True) # lowpass / downsample\n", " IQ[j,k,:] = sig.sosfilt(lpsos, x) # more low pass" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(3,1, sharex=True, figsize=(15,20))\n", "tt = np.linspace(0, b/(f_s), int(b/down_fac))\n", "ax = ax.flatten()\n", "\n", "p = np.where(tt/60 > 5)\n", "\n", "for m in range(len(freqs)):\n", " for n in range(len(y)):\n", " \n", " l, = ax[m].plot(tt[p]/60, np.real(IQ[n,m,p][0]), alpha=0.75, label=chans[chan_idx[n]] + '_I')\n", " ax[m].plot(tt[p]/60, np.imag(IQ[n,m,p][0]), alpha=0.75, label=chans[chan_idx[n]] + '_Q', c=l.get_color(), ls='--')\n", " ax[m].legend(ncol=2)\n", " ax[m].set_title(exc[m] + ' excitation = ' + str(freqs[m]) + ' Hz')\n", "\n", "plt.xlabel('Time [m]')\n", "plt.suptitle('UTC start time = ' + times)\n", "plt.savefig('Data/' + 'drmiLSCmatrix.pdf', bbox_inches='tight')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Some references\n", "* https://www.lsc-group.phys.uwm.edu/daswg/projects/nds-client/doc/manual/ch04.html\n", "* http://wiki.scipy.org/NumPy_for_Matlab_Users\n", "* https://www.pythonforthelab.com/\n", "* http://nbviewer.ipython.org/github/mantaraya36/201A-ipython/tree/master/Audio%20Filters.ipynb\n", "* http://stackoverflow.com/questions/15961979/how-do-i-plot-a-spectrogram-the-same-way-that-pylabs-specgram-does\n", "* http://matplotlib.org/api/pyplot_api.html#matplotlib.pyplot.imshow" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "a,b = y.shape\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ax.flatten()" ] }, { "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.7.1" } }, "nbformat": 4, "nbformat_minor": 2 }