sonification.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "749b0ddf",
   "metadata": {},
   "source": [
    "<img align=\"left\" src = https://project.lsst.org/sites/default/files/Rubin-O-Logo_0.png width=170 style=\"padding: 10px\"> \n",
    "<br><b>Little Demo: Data Sonification</b> <br>\n",
    "Contact author: Andrés A. Plazas Malagón <br>\n",
    "Last verified to run: 2024-07-31 <br>\n",
    "LSST Science Pipelines version: Weekly 2024_16 <br>\n",
    "Container size: small <be>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3f776002-eae8-426a-af45-d99988b06a65",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "This notebook demonstrates the use of the software [`Astronify`](https://astronify.readthedocs.io/en/latest/astronify/index.html) on the Rubin Science Platform (RSP) for sonifying astronomical data. `Astronify` is a toolkit developed by members of the Space Telescope Science Institute and designed to convert data series into sound by mapping one column to time and another to pitch. Usually, the observation time is scaled to listening time and flux is mapped to pitch, which is related to the sound frequency. Other software for sonification in astronomy, such as the one developed by the [`Rubin Rhapsodies` project](https://github.com/RileyWClarke/RubinRhapsodies/tree/main), map additional features such as flux error to volume (related to amplitude) and filter band to timbre.\n",
    "\n",
    "The `Astronify` package is under active development, with plans to expand its sonification capabilities further.\n",
    "\n",
    "Currently, there are issues with installing `Astronify` directly on the RSP due to dependencies like [`pyo`](https://github.com/belangeo/pyo) and [`libsndfile`](https://github.com/libsndfile/libsndfile). However, it is possible to work around this limitation by directly copying the necessary classes and functions from `Astronify`. \n",
    "\n",
    "This demo will feature a Kepler light curve, as demonstrated in the `Astronify` documentation, and a light curve from the Data Preview 0.2 (DP0.2) data, as shown in the `07a_DiaObject_samples` tutorial."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3b5324aa-5790-4a15-9ce3-96ce79414ac2",
   "metadata": {},
   "source": [
    "### 1.1 Import Packages\n",
    "\n",
    "The [`lightkurve`](https://docs.lightkurve.org/about/install.html) package can be locally installed by opening a terminal in the RSP and typing:\n",
    "\n",
    "- `python3 -m pip install lightkurve`\n",
    "\n",
    "The `astropy`, `warnings`, and `inspect` imports are required by the `Astronify` classes and functions used in this notebook and defined below in section 1.2. \n",
    "\n",
    "The `numpy` and `wave` imports are used to replace the `write()` method of `Astronify`'s `SoniSeries` class, which originally uses `pyo`, and write a `WAV` file with the audio to disk.  \n",
    "\n",
    "The `IPython` package is used to read and reproduce the `WAV` audio file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "57ca9517-8720-45c0-9b1f-6f737e8a9cfd",
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.utils.exceptions import AstropyWarning\n",
    "import lightkurve\n",
    "from astropy.table import Table, MaskedColumn\n",
    "from astropy.time import Time\n",
    "from astropy.visualization import (SqrtStretch, LogStretch, AsinhStretch, SinhStretch,\n",
    "                                   LinearStretch, MinMaxInterval, ManualInterval,\n",
    "                                   AsymmetricPercentileInterval)\n",
    "import numpy as np\n",
    "import warnings\n",
    "from inspect import signature, Parameter\n",
    "import wave\n",
    "import IPython\n",
    "\n",
    "from lsst.rsp import get_tap_service"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fdf314a6-7e0f-48b9-be6f-cd43cc89c36b",
   "metadata": {},
   "source": [
    "### 1.2 Define functions and parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "866cb1d4-7704-474a-ace4-9ab574cf1004",
   "metadata": {},
   "source": [
    "Import utilities from `Astronify` files [exceptions.py](https://github.com/spacetelescope/astronify/blob/main/astronify/utils/exceptions.py) and [pitch_mapping.py](https://github.com/spacetelescope/astronify/blob/main/astronify/utils/pitch_mapping.py) from the `Astronify` [repository](https://github.com/spacetelescope/astronify). Click on each cell to see the hidden code in the following two cells."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9abd30f1-67a8-4a8e-97dc-f528e6ef833b",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "class InvalidInputError(Exception):\n",
    "    \"\"\"\n",
    "    Exception to be issued when user input is incorrect in a \n",
    "    way that prevents the function from running.\n",
    "    \"\"\"\n",
    "    pass\n",
    "\n",
    "\n",
    "class InputWarning(AstropyWarning):\n",
    "    \"\"\"\n",
    "    Warning to be issued when user input is incorrect in\n",
    "    some way but doesn't prevent the function from running.\n",
    "    \"\"\"\n",
    "    pass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7a6b5399-26ae-4e1c-b1d2-d6257fc3165d",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "def data_to_pitch(data_array, pitch_range=[100, 10000], center_pitch=440, zero_point=\"median\",\n",
    "                  stretch='linear', minmax_percent=None, minmax_value=None, invert=False):\n",
    "    \"\"\"\n",
    "    Map data array to audible pitches in the given range, and apply stretch and scaling\n",
    "    as required.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    data_array : array-like\n",
    "        Data to map to pitch values. Individual data values should be floats.\n",
    "    pitch_range : array\n",
    "        Optional, default [100,10000]. Range of acceptable pitches in Hz.\n",
    "    center_pitch : float\n",
    "        Optional, default 440. The pitch in Hz where that the the zero point of the\n",
    "        data will be mapped to.\n",
    "    zero_point : str or float\n",
    "        Optional, default \"median\". The data value that will be mapped to the center\n",
    "        pitch. Options are mean, median, or a specified data value (float).\n",
    "    stretch : str\n",
    "        Optional, default 'linear'. The stretch to apply to the data array.\n",
    "        Valid values are: asinh, sinh, sqrt, log, linear\n",
    "    minmax_percent : array\n",
    "        Optional. Interval based on a keeping a specified fraction of data values\n",
    "        (can be asymmetric) when scaling the data. The format is [lower percentile,\n",
    "        upper percentile], where data values below the lower percentile and above the upper\n",
    "        percentile are clipped. Only one of minmax_percent and minmax_value should be specified.\n",
    "    minmax_value : array\n",
    "        Optional. Interval based on user-specified data values when scaling the data array.\n",
    "        The format is [min value, max value], where data values below the min value and above\n",
    "        the max value are clipped.\n",
    "        Only one of minmax_percent and minmax_value should be specified.\n",
    "    invert : bool\n",
    "        Optional, default False.  If True the pitch array is inverted (low pitches become high\n",
    "        and vice versa).\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    response : array\n",
    "        The normalized data array, with values in given pitch range.\n",
    "    \"\"\"\n",
    "    # Parsing the zero point\n",
    "    if zero_point in (\"med\", \"median\"):\n",
    "        zero_point = np.median(data_array)\n",
    "    if zero_point in (\"ave\", \"mean\", \"average\"):\n",
    "        zero_point = np.mean(data_array)\n",
    "\n",
    "    # The center pitch cannot be >= max() pitch range, or <= min() of pitch range.\n",
    "    # If it is, fall back to using the mean of the pitch range provided.\n",
    "    if center_pitch <= pitch_range[0] or center_pitch >= pitch_range[1]:\n",
    "        warnings.warn(\"Given center pitch is outside the pitch range, defaulting to the mean.\",\n",
    "                      InputWarning)\n",
    "        center_pitch = np.mean(pitch_range)\n",
    "\n",
    "    if (data_array == zero_point).all():  # All values are the same, no more calculation needed\n",
    "        return np.full(len(data_array), center_pitch)\n",
    "\n",
    "    # Normalizing the data_array and adding the zero point (so it can go through the same transform)\n",
    "    data_array = np.append(np.array(data_array), zero_point)\n",
    "\n",
    "    # Setting up the transform with the stretch\n",
    "    if stretch == 'asinh':\n",
    "        transform = AsinhStretch()\n",
    "    elif stretch == 'sinh':\n",
    "        transform = SinhStretch()\n",
    "    elif stretch == 'sqrt':\n",
    "        transform = SqrtStretch()\n",
    "    elif stretch == 'log':\n",
    "        transform = LogStretch()\n",
    "    elif stretch == 'linear':\n",
    "        transform = LinearStretch()\n",
    "    else:\n",
    "        raise InvalidInputError(\"Stretch {} is not supported!\".format(stretch))\n",
    "\n",
    "    # Adding the scaling to the transform\n",
    "    if minmax_percent is not None:\n",
    "        transform += AsymmetricPercentileInterval(*minmax_percent)\n",
    "\n",
    "        if minmax_value is not None:\n",
    "            warnings.warn(\"Both minmax_percent and minmax_value are set, minmax_value will be ignored.\",\n",
    "                          InputWarning)\n",
    "    elif minmax_value is not None:\n",
    "        transform += ManualInterval(*minmax_value)\n",
    "    else:  # Default, scale the entire image range to [0,1]\n",
    "        transform += MinMaxInterval()\n",
    "\n",
    "    # Performing the transform and then putting it into the pich range\n",
    "    pitch_array = transform(data_array)\n",
    "\n",
    "    if invert:\n",
    "        pitch_array = 1 - pitch_array\n",
    "\n",
    "    zero_point = pitch_array[-1]\n",
    "    pitch_array = pitch_array[:-1]\n",
    "\n",
    "    # In rare cases, the zero-point at this stage might be 0.0.\n",
    "    # One example is an input array of two values where the median() is the same as the\n",
    "    # lowest of the two values. In this case, the zero-point is 0.0 and will lead to error\n",
    "    # (divide by zero). Change to small value to avoid dividing by zero (in reality the choice\n",
    "    # of zero-point calculation by the user was probably poor, but not in purview to mandate or\n",
    "    # change user's choice here.  May want to consider providing info back to the user about the\n",
    "    # distribution of pitches actually used based on their sonification options in some way.\n",
    "    if zero_point == 0.0:\n",
    "        zero_point = 1E-6\n",
    "\n",
    "    if ((1/zero_point)*(center_pitch - pitch_range[0]) + pitch_range[0]) <= pitch_range[1]:\n",
    "        pitch_array = (pitch_array/zero_point)*(center_pitch - pitch_range[0]) + pitch_range[0]\n",
    "    else:\n",
    "        pitch_array = (((pitch_array-zero_point)/(1-zero_point))*(pitch_range[1] - center_pitch) +\n",
    "                       center_pitch)\n",
    "\n",
    "    return pitch_array"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c31ee49-dae1-4dd0-9075-15b45e460bee",
   "metadata": {},
   "source": [
    "Import classes from `Astronify` file [series.py](https://github.com/spacetelescope/astronify/blob/main/astronify/series/series.py) from the `Astronify` [repository](https://github.com/spacetelescope/astronify). The `write` method in the `SoniSeries` class has been replaced with a version that utilizes `numpy` and `wave` instead of `pyo`. Additionally, the `self._init_pyo()` method in the `SoniSeries` constructor is commented out to avoid dependency issues. Click on each cell to see the hidden code in the following two cells."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dccfef87-8770-4d14-88a4-c9f9ee21ea34",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "class PitchMap():\n",
    "\n",
    "    def __init__(self, pitch_func=data_to_pitch, **pitch_args):\n",
    "        \"\"\"\n",
    "        Class that encapsulates the data value to pitch function \n",
    "        and associated arguments.\n",
    "\n",
    "        Parameters\n",
    "        ----------\n",
    "        pitch_func : function\n",
    "            Optional. Defaults to `~astronify.utils.data_to_pitch`.\n",
    "            If supplying a function it should take a data array as the first\n",
    "            parameter, and all other parameters should be optional.\n",
    "        **pitch_args \n",
    "            Default parameters and values for the pitch function. Should include\n",
    "            all necessary arguments other than the data values.\n",
    "        \"\"\"\n",
    "\n",
    "        # Setting up the default arguments\n",
    "        if (not pitch_args) and (pitch_func == data_to_pitch):\n",
    "            pitch_args = {\"pitch_range\": [100, 10000],\n",
    "                          \"center_pitch\": 440,\n",
    "                          \"zero_point\": \"median\",\n",
    "                          \"stretch\": \"linear\"}\n",
    "        \n",
    "        self.pitch_map_func = pitch_func\n",
    "        self.pitch_map_args = pitch_args\n",
    "\n",
    "        \n",
    "    def _check_func_args(self):\n",
    "        \"\"\"\n",
    "        Make sure the pitch mapping function and argument dictionary match.\n",
    "\n",
    "        Note: This function does not check the the function gets all the required arguments.\n",
    "        \"\"\"\n",
    "        # Only test if both pitch func and args are set\n",
    "        if hasattr(self, \"pitch_map_func\") and hasattr(self, \"pitch_map_args\"):\n",
    "\n",
    "            # Only check parameters if there is no kwargs argument\n",
    "            param_types = [x.kind for x in signature(self.pitch_map_func).parameters.values()]\n",
    "            if Parameter.VAR_KEYWORD not in param_types:\n",
    "                for arg_name in list(self.pitch_map_args):\n",
    "                    if arg_name not in signature(self.pitch_map_func).parameters:\n",
    "                        wstr = \"{} is not accepted by the pitch mapping function and will be ignored\".format(arg_name)\n",
    "                        warnings.warn(wstr, InputWarning)\n",
    "                        del self.pitch_map_args[arg_name]\n",
    "\n",
    "    def __call__(self, data):\n",
    "        \"\"\"\n",
    "        Where does this show up?\n",
    "        \"\"\"\n",
    "        self._check_func_args()\n",
    "        return self.pitch_map_func(data, **self.pitch_map_args)\n",
    "\n",
    "    @property\n",
    "    def pitch_map_func(self):\n",
    "        \"\"\"\n",
    "        The pitch mapping function. \n",
    "        \"\"\"\n",
    "        return self._pitch_map_func\n",
    "\n",
    "    @pitch_map_func.setter\n",
    "    def pitch_map_func(self, new_func):\n",
    "        assert callable(new_func), \"Pitch mapping function must be a function.\"\n",
    "        self._pitch_map_func = new_func\n",
    "        self._check_func_args()\n",
    "\n",
    "    @property\n",
    "    def pitch_map_args(self):\n",
    "        \"\"\"\n",
    "        Dictionary of additional arguments (other than the data array)\n",
    "        for the pitch mapping function.\n",
    "        \"\"\"\n",
    "        return self._pitch_map_args\n",
    "\n",
    "    @pitch_map_args.setter\n",
    "    def pitch_map_args(self, new_args):\n",
    "        assert isinstance(new_args, dict), \"Pitch mapping function args must be in a dictionary.\"\n",
    "        self._pitch_map_args = new_args\n",
    "        self._check_func_args()\n",
    "\n",
    "             \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "59e100fb-1156-4159-84c5-9553fee794cf",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "class SoniSeries():\n",
    "\n",
    "    def __init__(self, data, time_col=\"time\", val_col=\"flux\"):\n",
    "        \"\"\"\n",
    "        Class that encapsulates a sonified data series.\n",
    "\n",
    "        Parameters\n",
    "        ----------\n",
    "        data : `astropy.table.Table`\n",
    "            The table of data to be sonified.\n",
    "        time_col : str\n",
    "            Optional, default \"time\". The data column to be mapped to time.\n",
    "        val_col : str\n",
    "            Optional, default \"flux\". The data column to be mapped to pitch.\n",
    "        \"\"\"\n",
    "        self.time_col = time_col\n",
    "        self.val_col = val_col\n",
    "        self.data = data\n",
    "\n",
    "        # Default specs\n",
    "        self.note_duration = 0.5  # note duration in seconds\n",
    "        self.note_spacing = 0.01  # spacing between notes in seconds\n",
    "        self.gain = 0.05  # default gain in the generated sine wave. pyo multiplier, -1 to 1.\n",
    "        self.pitch_mapper = PitchMap(data_to_pitch)\n",
    "\n",
    "        # self._init_pyo()\n",
    "\n",
    "    def _init_pyo(self):\n",
    "        self.server = pyo.Server()\n",
    "        self.streams = None\n",
    "\n",
    "    @property\n",
    "    def data(self):\n",
    "        \"\"\" The data table (~astropy.table.Table). \"\"\"\n",
    "        return self._data\n",
    "\n",
    "    @data.setter\n",
    "    def data(self, data_table):\n",
    "\n",
    "        if not isinstance(data_table, Table):\n",
    "            raise TypeError('Data must be an astropy.table.Table object.')\n",
    "\n",
    "        for c in list(data_table.columns):\n",
    "            data_table.rename_column(c, c.lower())\n",
    "\n",
    "\n",
    "        if self.time_col not in data_table.columns:\n",
    "            raise AttributeError(f\"Input Table must contain time column '{self.time_col}'\")\n",
    "\n",
    "        if self.val_col not in data_table.columns:\n",
    "            raise AttributeError(f\"Input Table must contain a value column '{self.val_col}'\")\n",
    "\n",
    "        # Removing any masked values as they interfere with the sonification\n",
    "        if isinstance(data_table[self.val_col], MaskedColumn):\n",
    "            data_table = data_table[~data_table[self.val_col].mask]\n",
    "        if isinstance(data_table[self.time_col], MaskedColumn):\n",
    "            data_table = data_table[~data_table[self.time_col].mask]\n",
    "\n",
    "        # Removing any nans as they interfere with the sonification\n",
    "        data_table = data_table[~np.isnan(data_table[self.val_col])]\n",
    "\n",
    "        # making sure we have a float column for time\n",
    "        if isinstance(data_table[self.time_col], Time):\n",
    "            float_col = \"asf_time\"\n",
    "            data_table[float_col] = data_table[self.time_col].jd\n",
    "            self.time_col = float_col\n",
    "            \n",
    "        self._data = data_table\n",
    "\n",
    "    @property\n",
    "    def time_col(self):\n",
    "        \"\"\" The data column mappend to time when sonifying. \"\"\"\n",
    "        return self._time_col\n",
    "\n",
    "    @time_col.setter\n",
    "    def time_col(self, value):\n",
    "        assert isinstance(value, str), 'Time column name must be a string.'\n",
    "        self._time_col = value\n",
    "\n",
    "    @property\n",
    "    def val_col(self):\n",
    "        \"\"\" The data column mappend to putch when sonifying. \"\"\"\n",
    "        return self._val_col\n",
    "\n",
    "    @val_col.setter\n",
    "    def val_col(self, value):\n",
    "        assert isinstance(value, str), 'Value column name must be a string.'\n",
    "        self._val_col = value\n",
    "\n",
    "    @property\n",
    "    def pitch_mapper(self):\n",
    "        \"\"\" The pitch mapping object that takes data values to pitch values (Hz). \"\"\"\n",
    "        return self._pitch_mapper\n",
    "\n",
    "    @pitch_mapper.setter\n",
    "    def pitch_mapper(self, value):\n",
    "        self._pitch_mapper = value\n",
    "\n",
    "    @property\n",
    "    def gain(self):\n",
    "        \"\"\" Adjustable gain for output. \"\"\"\n",
    "        return self._gain\n",
    "\n",
    "    @gain.setter\n",
    "    def gain(self, value):\n",
    "        self._gain = value\n",
    "\n",
    "    @property\n",
    "    def note_duration(self):\n",
    "        \"\"\" How long each individual note will be in seconds.\"\"\"\n",
    "        return self._note_duration\n",
    "\n",
    "    @note_duration.setter\n",
    "    def note_duration(self, value):\n",
    "        # Add in min value check\n",
    "        self._note_duration = value\n",
    "\n",
    "    @property\n",
    "    def note_spacing(self):\n",
    "        \"\"\" The spacing of the notes on average (will adjust based on time) in seconds. \"\"\"\n",
    "        return self._note_spacing\n",
    "\n",
    "    @note_spacing.setter\n",
    "    def note_spacing(self, value):\n",
    "        # Add in min value check\n",
    "        self._note_spacing = value\n",
    "        \n",
    "    def sonify(self):\n",
    "        \"\"\"\n",
    "        Perform the sonification, two columns will be added to the data table: asf_pitch, and asf_onsets. \n",
    "        The asf_pitch column will contain the sonified data in Hz.\n",
    "        The asf_onsets column will contain the start time for each note in seconds from the first note.\n",
    "        Metadata will also be added to the table giving information about the duration and spacing \n",
    "        of the sonified pitches, as well as an adjustable gain.\n",
    "        \"\"\"\n",
    "        data = self.data\n",
    "        exptime = np.median(np.diff(data[self.time_col]))\n",
    "\n",
    "        data.meta[\"asf_exposure_time\"] = exptime\n",
    "        data.meta[\"asf_note_duration\"] = self.note_duration\n",
    "        data.meta[\"asf_spacing\"] = self.note_spacing\n",
    "        \n",
    "        data[\"asf_pitch\"] = self.pitch_mapper(data[self.val_col])\n",
    "        data[\"asf_onsets\"] = [x for x in (data[self.time_col] - data[self.time_col][0])/exptime*self.note_spacing]\n",
    "\n",
    "    def play(self):\n",
    "        \"\"\"\n",
    "        Play the data sonification.\n",
    "        \"\"\"\n",
    "\n",
    "        # Making sure we have a clean server\n",
    "        if self.server.getIsBooted():\n",
    "            self.server.shutdown()\n",
    "\n",
    "        self.server.boot()\n",
    "        self.server.start()\n",
    "\n",
    "        # Getting data ready\n",
    "        duration = self.data.meta[\"asf_note_duration\"]\n",
    "        pitches = np.repeat(self.data[\"asf_pitch\"], 2)\n",
    "        delays = np.repeat(self.data[\"asf_onsets\"], 2)\n",
    "\n",
    "        # TODO: This doesn't seem like the best way to do this, but I don't know\n",
    "        # how to make it better\n",
    "        env = pyo.Linseg(list=[(0, 0), (0.01, 1), (duration - 0.1, 1),\n",
    "                               (duration - 0.05, 0.5), (duration - 0.005, 0)],\n",
    "                         mul=[self.gain for i in range(len(pitches))]).play(\n",
    "                             delay=list(delays), dur=duration)\n",
    "\n",
    "        self.streams = pyo.Sine(list(pitches), 0, env).out(delay=list(delays),\n",
    "                                                           dur=duration)\n",
    "\n",
    "    def stop(self):\n",
    "        \"\"\"\n",
    "        Stop playing the data sonification.\n",
    "        \"\"\"\n",
    "        self.streams.stop() \n",
    "\n",
    "    def write(self, filepath):\n",
    "        \"\"\"\n",
    "        Save data sonification to the given file in WAV format.\n",
    "    \n",
    "        Parameters\n",
    "        ----------\n",
    "        filepath : str\n",
    "            The path to the output file.\n",
    "        \"\"\"\n",
    "        # Sampling parameters\n",
    "        sample_rate = 44100  # Standard sample rate in Hz\n",
    "        max_amplitude = 32767  # Max amplitude for 16-bit audio\n",
    "        amplitude = 0.5  # Amplitude scaling factor\n",
    "    \n",
    "        # Prepare the data for writing\n",
    "        duration = self.data.meta[\"asf_note_duration\"]\n",
    "        pitches = self.data[\"asf_pitch\"]\n",
    "        onsets = self.data[\"asf_onsets\"]\n",
    "\n",
    "        # Calculate the total number of samples\n",
    "        total_duration = onsets[-1] + duration\n",
    "        total_samples = int(sample_rate * total_duration)\n",
    "    \n",
    "        # Create an array to hold the audio data\n",
    "        audio_data = np.zeros(total_samples)\n",
    "    \n",
    "        for pitch, onset in zip(pitches, onsets):\n",
    "            if pitch <= 0:\n",
    "                continue  # Skip silence or invalid pitch values\n",
    "    \n",
    "            start_sample = int(sample_rate * onset)\n",
    "            end_sample = start_sample + int(sample_rate * duration)\n",
    "            t = np.linspace(0, duration, end_sample - start_sample, endpoint=False)\n",
    "            sine_wave = amplitude * np.sin(2 * np.pi * pitch * t)\n",
    "    \n",
    "            # Apply a simple fade-in and fade-out envelope\n",
    "            fade_length = int(sample_rate * 0.01)  # 10 ms fade\n",
    "            envelope = np.ones_like(sine_wave)\n",
    "            envelope[:fade_length] = np.linspace(0, 1, fade_length)\n",
    "            envelope[-fade_length:] = np.linspace(1, 0, fade_length)\n",
    "            sine_wave *= envelope\n",
    "            \n",
    "            # Ensure the generated samples fit within the audio_data array\n",
    "            audio_data[start_sample:end_sample] += sine_wave\n",
    "    \n",
    "        # Normalize the audio data to prevent clipping\n",
    "        audio_data = (audio_data / np.max(np.abs(audio_data))) * max_amplitude\n",
    "        audio_data = audio_data.astype(np.int16)  # Convert to 16-bit PCM format\n",
    "    \n",
    "        # Write to WAV file\n",
    "        with wave.open(filepath, 'w') as wf:\n",
    "            wf.setnchannels(1)  # Mono sound\n",
    "            wf.setsampwidth(2)  # 2 bytes per sample (16 bits)\n",
    "            wf.setframerate(sample_rate)\n",
    "            wf.writeframes(audio_data.tobytes())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8a02c273-7d1d-4b6c-95e7-9cdeec318ef2",
   "metadata": {},
   "source": [
    "## 2. Sonification of a Kepler lightcurve\n",
    "\n",
    "Reproduce one example from the `Astronify` documentation."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b059874d-cb2c-42bf-86d4-21e499a63089",
   "metadata": {},
   "source": [
    "The following query is modified slightly from what is shown in the [documentation](https://astronify.readthedocs.io/en/latest/astronify/index.html) to avoid deprecation messages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f7b22d35-5797-441c-ad88-994e727b0dd1",
   "metadata": {},
   "outputs": [],
   "source": [
    "kep12b_lc = lightkurve.search_lightcurve(\"KIC 11804465\", cadence=\"long\", quarter=1).download_all()[0].to_table()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb8b0638-47ec-46aa-a346-b0791855e685",
   "metadata": {},
   "source": [
    "Print the table to see its format. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "01201fa4-e457-4897-a8d5-c1dd3a6e98d7",
   "metadata": {},
   "outputs": [],
   "source": [
    "print (kep12b_lc)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bb1b7d68-32b9-4f6b-8f11-c1ccba8baca9",
   "metadata": {},
   "source": [
    "By default, `Soniseries` looks for columns labeled \"time\" and \"flux\", but this can be changed in the argument of the class (see the [documentation](https://astronify.readthedocs.io/en/latest/astronify/index.html))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cb13b391-fa79-490b-82cf-89620e4e19c1",
   "metadata": {},
   "outputs": [],
   "source": [
    "kep12b_obj = SoniSeries(kep12b_lc)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1ab833f7-b6f3-4a0f-9ae8-4fee166094f8",
   "metadata": {},
   "source": [
    "Sonify the data. See the `Astronify` [documentation](https://astronify.readthedocs.io/en/latest/astronify/index.html) for more options and a description of the algorithm used.\n",
    "\n",
    "The warning can be ignored."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0368283c-1d86-4473-af23-5733e6750cad",
   "metadata": {},
   "outputs": [],
   "source": [
    "kep12b_obj.sonify() "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5f05402d-e75b-4679-a52b-ee22c59b5924",
   "metadata": {},
   "source": [
    "The `SoniSeries` class has a `play()` method that is used to play the sonified data. This method cannot be used on the RSP, as the platform does not support direct sound playback.\n",
    "\n",
    "Instead, the sonified data can be written to a file in the RSP (such as a `WAV` file) and played in the notebook using `IPython`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a4c34603-8baa-41bc-85cc-a91f807b4ffc",
   "metadata": {},
   "source": [
    "Use the redifined `write` function to produce an audiofile in `WAV` format."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8a04a349-c1c2-40b0-8cf8-d993ed39b30b",
   "metadata": {},
   "outputs": [],
   "source": [
    "kep12b_obj.write(\"./kepler12_sonification_demo.wav\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4921b399-95d4-420b-a0ff-fd0821ba37b4",
   "metadata": {},
   "source": [
    "Play the file with `IPython`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "53cba4d9-9ce7-402e-8d04-8114943da41e",
   "metadata": {},
   "outputs": [],
   "source": [
    "IPython.display.Audio(filename=\"./kepler12_sonification_demo.wav\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ac57190-c550-45be-b41c-02bce9e34616",
   "metadata": {},
   "source": [
    "## 3. Sonification of a DP0.2 light curve.\n",
    "\n",
    "Use data from DP0.2 tutorial `07a_DiaObject_samples`, section 4.2.2 \"Light curve for epochs near SNR>5 detections\", to sonify the light curve of `DiaObject` with ID `1250953961339360185`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2c813697-c585-451d-b710-8016de757249",
   "metadata": {},
   "outputs": [],
   "source": [
    "DiaObjID = 1250953961339360185\n",
    "print(DiaObjID)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9a2abf8d-e30d-45ad-862c-a07534413d1d",
   "metadata": {},
   "source": [
    "Use the time and difference-image point source (PS) flux data from tutorial `07a_DiaObject_samples`. Use the `DiaSrc['midPointTai']` and `DiaSrc['psFlux']` data in the `r` band from section 4.2.2 of that tutorial.  Time is reported in the DiaSource table as `midPointTai`, which is in the SI unit of \"TAI\" (<a href=\"https://en.wikipedia.org/wiki/International_Atomic_Time\">International Atomic Time</a>), and is presented in days (in particular, as \"<a href=\"https://en.wikipedia.org/wiki/Julian_day\">Modified Julian Days</a>\").\n",
    "\n",
    "The query is a simplified version of the query in that section of the tutorial, retrieving only the quantities needed for this demo."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0462fd00-b05f-46c0-b86b-1429e78f44b4",
   "metadata": {},
   "outputs": [],
   "source": [
    "service = get_tap_service(\"tap\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "29121b11-3b60-4a7d-8fdc-239ea448872c",
   "metadata": {},
   "outputs": [],
   "source": [
    "results = service.search(\"SELECT psFlux, \"\n",
    "                         \"filterName, midPointTai \"\n",
    "                         \"FROM dp02_dc2_catalogs.DiaSource \"\n",
    "                         \"WHERE diaObjectId = \"+str(DiaObjID))\n",
    "DiaSrc = results.to_table()\n",
    "del results"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4675b408-2b13-4d3a-9d32-9b9b1db1bcc6",
   "metadata": {},
   "source": [
    "Get the time and flux in the `r` band."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "48b279e0-d493-45d6-adb2-1985a6ce3d5a",
   "metadata": {},
   "outputs": [],
   "source": [
    "filt = 'r'\n",
    "fx = np.where(DiaSrc['filterName'] == filt)[0]\n",
    "time_, flux_ = DiaSrc['midPointTai'][fx], DiaSrc['psFlux'][fx]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5010e46a-e7b3-412b-a4b6-748d5db9eee3",
   "metadata": {},
   "source": [
    "`SoniSeries` expect the data to be sorted by time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "33b0cdc5-de94-446e-beb9-0af56970acbf",
   "metadata": {},
   "outputs": [],
   "source": [
    "index = np.argsort(time_)\n",
    "time = time_[index]\n",
    "flux = flux_[index]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "45a42a11-c545-4004-b4fa-722bc3a0307c",
   "metadata": {},
   "source": [
    "Create an `astropy` Table to instantiate the `SoniSeries` object, and call the `sonify` mehod. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b9d91d73-d866-4219-933c-cfab1e8bcec7",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_table = Table({\"time\": time,\n",
    "                    \"flux\": flux})\n",
    "\n",
    "data_soni = SoniSeries(data_table)\n",
    "data_soni.sonify()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "828ee39b-4918-4c3b-a8fa-317d90b09e00",
   "metadata": {},
   "source": [
    "Write the data into a `WAV` file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "705719c7-c2b3-426f-9790-6f86b15bbe9b",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_soni.write(\"./dia_object_sonification_demo.wav\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5354550d-93e2-46a4-bf8a-9a0c2d651527",
   "metadata": {},
   "source": [
    "Reproduce the audio file. The sounds are sparsers compared to the more constant sounds from the Kepler light curve from section 2, which has a higher cadence. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "087a74cb-ed87-4f49-a7ba-3dc391ffd3d7",
   "metadata": {},
   "outputs": [],
   "source": [
    "IPython.display.Audio(filename=\"./dia_object_sonification_demo.wav\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e19cbf0d-916b-4f46-a7cd-09c05e89ee6b",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "LSST",
   "language": "python",
   "name": "lsst"
  },
  "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.11.7"
  },
  "toc-autonumbering": false
 },
 "nbformat": 4,
 "nbformat_minor": 5
}