Match_Data_from_File_1vs1.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "#Codes by Shucheng Yang\n",
    "import numpy as np\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import pycbc.psd\n",
    "\n",
    "from pycbc.types.timeseries import TimeSeries\n",
    "from pycbc.types.frequencyseries import FrequencySeries\n",
    "from pycbc import filter\n",
    "\n",
    "\n",
    "from scipy.interpolate import interp1d #导入scipy库\n",
    "import os\n",
    "\n",
    "\n",
    "def file2series(path):\n",
    "    with open(path, 'r') as f:\n",
    "        datafile = f.readlines()\n",
    "\n",
    "    datafile = [np.fromstring(item.replace(\"D\",\"e\").replace('\\n','').strip(), dtype= float, sep =' ') for item in datafile]\n",
    "    datafile = np.array(datafile)\n",
    "    ans = {}\n",
    "    \n",
    "    ans['timeVec'] = datafile[:,0]\n",
    "    ans['n'] = len(ans['timeVec'])\n",
    "    ans['sampIntrvl'] = ans['timeVec'][1]-ans['timeVec'][0]\n",
    "    \n",
    "    ans['hp'] = TimeSeries( datafile[:,5], delta_t=ans['sampIntrvl'], dtype = float,copy=True)\n",
    "    ans['hc'] = TimeSeries( datafile[:,6], delta_t=ans['sampIntrvl'], dtype = float,copy=True)    \n",
    "    return ans\n",
    "\n",
    "\n",
    "#绘制波形图\n",
    "def plotwave(hp,hc):\n",
    "    plt.figure(figsize=(8,8/1.5))\n",
    "    plt.plot(hp.sample_times, hp, label = '$h_{+}$')\n",
    "    plt.plot(hc.sample_times, hc, label = '$h_{\\\\times}$')\n",
    "    # plt.xlim(40,40.5)\n",
    "    # plt.ylim(- 2e-22,2e-22)\n",
    "    plt.xlabel(\"Time / s\")\n",
    "    plt.ylabel(\"$Strain$\")\n",
    "    plt.legend()\n",
    "    plt.tight_layout()\n",
    "    plt.grid(linestyle = \"dotted\", color = \"#d3d3d3\" , which=\"both\")\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "def stdf(h, f_max, freqIntrvl):\n",
    "    #interp1d\n",
    "    hf = h.to_frequencyseries()    \n",
    "    freqVec = hf.sample_frequencies\n",
    "    hf_value = np.array(hf)\n",
    "    fit=interp1d(freqVec, hf_value ,fill_value = \"extrapolate\")  \n",
    "    #\n",
    "    stdFreqVec = np.arange(0, f_max, freqIntrvl)\n",
    "    hf_value_std = fit(stdFreqVec)\n",
    "    stdf_series = pycbc.types.frequencyseries.FrequencySeries(hf_value_std, delta_f=freqIntrvl, epoch='', copy=True)\n",
    "    return stdf_series\n",
    "\n",
    "\n",
    "def overlap_match_func(hf1, hf2, psd, f_min, f_max):\n",
    "    roverlap = filter.matchedfilter.overlap(hf1, hf2, psd=psd, low_frequency_cutoff=f_min , high_frequency_cutoff=f_max, normalized=True)  \n",
    "    amplitude1 = filter.matchedfilter.sigmasq(hf1, psd=psd, low_frequency_cutoff=f_min, high_frequency_cutoff=f_max )\n",
    "    amplitude2 = filter.matchedfilter.sigmasq(hf2, psd=psd, low_frequency_cutoff=f_min, high_frequency_cutoff=f_max )\n",
    "    rmatch,nn =filter.matchedfilter.match(hf1, hf2, psd=psd, low_frequency_cutoff=f_min , high_frequency_cutoff=f_max,v1_norm=True ,v2_norm=True)\n",
    "    rmatch = rmatch / np.sqrt(amplitude1 * amplitude2)\n",
    "    return roverlap, rmatch"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "257022 0.000295529457\n",
      "212102 0.000295529457\n"
     ]
    }
   ],
   "source": [
    "path1 = \"/home/ysc/Desktop/nw1/heco60_NS14_a60s50_R9_nu250_f10Hz_fisco.dat\"\n",
    "ans1 = file2series(path1)\n",
    "\n",
    "hp1 = ans1[\"hp\"]\n",
    "hc1 = ans1[\"hc\"]\n",
    "timeVec1 = ans1[\"timeVec\"]\n",
    "n1 = len(hc1)\n",
    "sampIntrvl1 = ans1[\"sampIntrvl\"]\n",
    "print(n1,sampIntrvl1)\n",
    "\n",
    "path2 = \"/home/ysc/Desktop/nw2/hbh60_NS14_a60_nu250_f10Hz_fisco.dat\"\n",
    "ans2 = file2series(path2)\n",
    "\n",
    "hp2 = ans2[\"hp\"]\n",
    "hc2 = ans2[\"hc\"]\n",
    "timeVec2 = ans2[\"timeVec\"]\n",
    "n2 = ans2[\"n\"]\n",
    "sampIntrvl2 = ans2[\"sampIntrvl\"]\n",
    "print(n2,sampIntrvl1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 参数设置"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "采样频率为3383.757444Hz，信号持续时间62s, 时域信号采样212102 个点\n",
      "信号中可分析最大频率为1691.878722Hz\n",
      "\n",
      "采样周期，即时域分辨率为0.000296s\n",
      "信号频域的频率间隔，即频域分辨率为0.015953Hz\n"
     ]
    }
   ],
   "source": [
    "assert sampIntrvl1 == sampIntrvl2 \n",
    "\n",
    "sampIntrvl = sampIntrvl1\n",
    "n_min = min(n1,n2)          #采样点数(Sampling Number), 有时也称为信号长度(Length of Signal) 2^16为2的幂时，快速傅里叶变化效率最高\n",
    "                            #n =  duration * sampFreqint = (duration / sampIntrvl)\n",
    "sampFreq = 1/sampIntrvl  #采样频率(Sampling frequency)，单位时间样本点个数，应大于 2f（即Nyquist频率)\n",
    "\n",
    "duration = n_min/sampFreq #40        #信号持续时间(duration of signal) 2^16, 0.75d, 2^25 ,1.06yr\n",
    "\n",
    "sampIntrvl = 1.0 / sampFreq                   #采样周期(Sampling period)，隔多少时间取样一次，或步长\n",
    "freqIntrvl =  1 / (n_min * sampIntrvl)                     #傅里叶变换 频率分辨率(Frequency Interval) \n",
    "                                              # freqIntrvl = 1 / duration = 1 / (n * sampIntrvl)\n",
    "                                              #            = sampFreq / n  \n",
    "        \n",
    "f_min = 20              #低于此频率的psd将被设置为0\n",
    "f_max = sampFreq/2             #信号模式的最大频率\n",
    "\n",
    "print(\"采样频率为%fHz，信号持续时间%ds, 时域信号采样%d 个点\"%(sampFreq,duration,n_min))\n",
    "print(\"信号中可分析最大频率为%fHz\"%f_max)\n",
    "print(\"\\n采样周期，即时域分辨率为%fs\"%(sampIntrvl))\n",
    "print(\"信号频域的频率间隔，即频域分辨率为%fHz\"%freqIntrvl);\n",
    "\n",
    "\n",
    "# 列出lalsuite内置的解析psd (没发现有LISA的，下面以LIGO的作为示例)\n",
    "# print(pycbc.psd.get_lalsim_psd_list())\n",
    "\n",
    "#示例，psd参见， https://dcc.ligo.org/LIGO-T1800044/public\n",
    "psd = pycbc.psd.from_string('aLIGOaLIGODesignSensitivityT1800044', n_min , freqIntrvl, f_min)\n",
    "\n",
    "# #绘制 频率 - sqrt(PSD) 图\n",
    "# plt.figure(figsize=(8,8/1.5))\n",
    "# plt.loglog(psd.sample_frequencies, np.sqrt(psd), label = 'LIGO-psd')\n",
    "# plt.xlim(2e1,1e3)\n",
    "# plt.xlabel(\"Frequency / Hz\")\n",
    "# plt.ylabel(\"$\\sqrt{S_{n}(f) \\ / \\ (Hz^{-1})}$\")\n",
    "# plt.legend()\n",
    "# plt.tight_layout()\n",
    "# plt.grid(linestyle = \"dotted\", color = \"#d3d3d3\" , which=\"both\")\n",
    "# plt.show()\n",
    "\n",
    "# 波形图\n",
    "# plotwave(hp2,hc2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# match and overlap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Results:\n",
      "h+: match = 0.9939, overlap = 0.9704\n",
      "hx: match = 0.9922, overlap = 0.9705 \n",
      "\n"
     ]
    }
   ],
   "source": [
    "hpf1 = stdf(hp1, f_max, freqIntrvl)\n",
    "hpf2 = stdf(hp2, f_max, freqIntrvl)\n",
    "\n",
    "hcf1 = stdf(hc1, f_max, freqIntrvl)\n",
    "hcf2 = stdf(hc2, f_max, freqIntrvl)\n",
    "\n",
    "overlap_hp, match_hp = overlap_match_func(hpf1, hpf2, psd, f_min, f_max)\n",
    "overlap_hc, match_hc = overlap_match_func(hcf1, hcf2, psd, f_min, f_max)\n",
    "\n",
    "print(\"Results:\\nh+: match = {1:6.4f}, overlap = {0:6.4f}\".format(overlap_hp, match_hp))\n",
    "print(\"hx: match = {1:6.4f}, overlap = {0:6.4f} \\n\".format(overlap_hc, match_hc))"
   ]
  }
 ],
 "metadata": {
  "hide_input": false,
  "kernelspec": {
   "display_name": "pycbc",
   "language": "python",
   "name": "pycbc"
  },
  "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.3"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}