Sound Analysis-Static Figures.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Sound Analysis for Summer Ventures - Static Figures"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This version of the Sound Analysis notebook produces static figures using *Matplotlib*. These static figures are good for use in PowerPoint presentations and in papers. They can be exported to PNG, PDF, or JPEG files with the figure. There is also an *Interactive Figures* verion of this notebook."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This Jupyter notebook will allow you to do analysis of sound recordings and make graphs using the Plotly modules. Get started by following the directions before each cell. Change the parameters as needed, and execute the commands in the cell by entering `Shift-Enter` or `Shift-Return` with your the cell selected or with your cursor in the cell."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remember that you the order that you enter cells matters, not the order in which they appear in the notebook."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load Python Modules"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Enter the cell below to load the Python modules you will need. Add any additional modules you need for your specific calculations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np # Numerical analysis package\n",
    "import scipy.signal as signal # Signal processing module\n",
    "import matplotlib.pyplot as plt # Plotting module\n",
    "from SVSound import wavefile # Module to read WAV files"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read WAV File(s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You will use `wavefile.read` to read the information and data from sound recordings in the WAV format. This function was written to read WAV files written by many different devices, including the smartphone apps you may use."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `wavefile.read` function returns the file parameters and the data so you can use them for your analysis. As with most Python programs, enter `?wavefile` and `?wavefile.read` to get more information about them."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Enter the cell below to read the file `Beats.wav`, which has been loaded into this project. If you want to read a different WAV file, upload it to your project folder and change the `\"Beats.wav\"` in the command below to your filename (in quotes). (As always, feel free to change variable names `beats_info` and `beats_wave` to something that works for you.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "beats_info, beats_wave = wavefile.read(\"Beats.wav\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we are going to make variable with the times of each sample in `beats_wave`. We do this by taking a list of all the samples (0, 1, 2, 3, ...) amd dividing by the sampling frequency (because period $T = 1/f$). Do this by entering the cell below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "beats_times = np.arange(beats_info['Nsamples']) / beats_info['fs']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we will open a second WAV file, `OrganPipeOpenClosed.wav` so you can have two examples. Enter the cell below to run the same command we did above for this file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "pipe_info, pipe_wave = wavefile.read(\"OrganPipeOpenClosed.wav\")\n",
    "pipe_times = np.arange(pipe_info['Nsamples']) / pipe_info['fs']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot Waveforms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The cell below generates a plot of the entire `beats_wave` variable vs. the time values stored in `beats_times`. The resulting graph has buttons for zooming, panning, etc. Enter the cell below to beterate the graph."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=[12, 6])\n",
    "ax.plot(beats_times, beats_wave, '-b')\n",
    "ax.set_xlabel('Time (s)')\n",
    "ax.set_ylabel('Relative Pressure');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also save a static image of the plot for use in a report or presentation. You cannot zoom and pan a static image. Enter the cell below to save a PNG file containing the plot. You can also save a plot to JPEG, PDF, and SVG by changing the extension (part after the period) in the exported filename."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig.savefig(\"Beats Wave Plot.pdf\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now plot the second sound `pipe_wave`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=[12, 6])\n",
    "ax.plot(pipe_times, pipe_wave, '-b')\n",
    "ax.set_xlabel('Time (s)')\n",
    "ax.set_ylabel('Relative Pressure');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can save this plot using the same commands as the previous plot, replacing the name `beats_wavefig` with `pipe_wavfig`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Power Spectrum"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we will produce a power spectrum of the entire first sound stored in `beats_wave` with parameters stored in `beats_info`. Begin by displaying the sampling frequency."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is a tradeoff in FFTs and power spectra. You need to take more FFT samples to get a greater frequency resolution, but that requires you to sample a larger total time. If you are interested in frequency resolution, choose a large FFT size. If you want averaging, choose a smaller FFT size."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we will produce power spectra of portions of the first sound. First enter the number of FFT points you want to use. In this case, we are going to set the number of FFT points to the same number as the sampling frequency. This will give us a total sampling time of 1 s and a frequency resolution of 1 Hz."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "nfft = beats_info[\"fs\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now enter the cell below to calculate the power spectrum of the entire sound."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "beats_freq, beats_ps = signal.welch(beats_wave, beats_info['fs'], nperseg=nfft)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now enter the cell below to graph the power spectrum. Notice that we take the logarithm of the power spectrum and multiply by 10 before plotting. This converts it to decibels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=[12, 6])\n",
    "ax.plot(beats_freq, 10*np.log10(beats_ps), '-b')\n",
    "ax.set_xlabel('Frequency (Hz)')\n",
    "ax.set_ylabel('Power Spectral Density (dB)');\n",
    "ax.set_xlim([0, 2000])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we are going to find the power spectrum of the portion of the first sound before the beat with only one tuning fork playing. There are several ways to obtain a segment of the sound, including taking a \"slice\" of the `beats_wave` variable. The easiest way is just to read a segment of the soung because the wavefile.read function will pick out the times for you.\n",
    "\n",
    "Read the segment of Beats.wav from 1.65 s to 2.68 s by entering the following cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "onefork_info, onefork_wave = wavefile.read(\"Beats.wav\", t0=1.65, t1=2.68)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now take a power spectrum of this sound and plot it by entering the following cell. (Notice how we set the x-axis to only show frequencies from 0 Hz to 1000 Hz. You can still zoom out.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "nfft = onefork_info[\"fs\"]\n",
    "onefork_freq, onefork_ps = signal.welch(onefork_wave, onefork_info['fs'], nperseg=nfft)\n",
    "fig, ax = plt.subplots(figsize=[12, 6])\n",
    "ax.plot(onefork_freq, 10*np.log10(onefork_ps), '-b')\n",
    "ax.set_xlim([0, 1000])\n",
    "ax.set_ylim([-10, 80])\n",
    "ax.set_xlabel('Frequency (Hz)')\n",
    "ax.set_ylabel('Power Spectral Density (dB)');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we will take the power spectrum of the portion of the recording from 2.8 s to the end. Read this segment from the file by entering the following cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "twofork_info, twofork_wave = wavefile.read(\"Beats.wav\", t0=2.8, t1=-1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we will use the same parameters to make a power spectrum of this sound."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "nfft = twofork_info[\"fs\"]\n",
    "twofork_freq, twofork_ps = signal.welch(twofork_wave, twofork_info['fs'], nperseg=nfft)\n",
    "fig, ax = plt.subplots(figsize=[12, 6])\n",
    "ax.plot(twofork_freq, 10*np.log10(twofork_ps), '-b')\n",
    "ax.set_xlim([0, 1000])\n",
    "ax.set_ylim([-10, 80])\n",
    "ax.set_xlabel('Frequency (Hz)')\n",
    "ax.set_ylabel('Power Spectral Density (dB)');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Spectrogram"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another important tool for spectral analysis is the spectrogram. A spectrogram is a series of short-time power spectra, which shows how the frequency content changes in time. A spectrogram is a 3-D plot. One axis (usually the x-axis) is time. The second axis (usually the y-axis) is frequency. The z-axis (represented by color) is power spectral density, which represents how much signal power is in that frequency bin. (Note that sometimes but not often spectrograms have time on the y-axis and frequency on the x-axis.) Each slice of a spectrogram at a given time is a power spectrum."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Two important parameters for spectrograms are the number of FFT points. This is the number of samples used to calculate each power spectrum slice. Increasing the number of FFT samples increases the frequency resolution, but it decreases the time resolution because each power spectrum slice takes more time (more samples = more time). The second important parameter is the overlap between FFTs. This means the samples used for each FFT contain some of the samples from the previous FFT. The optimal overlap is half of the points used for each FFT.  In other words, it you use 2048 points for your FFT, the optimal overlap would be 1024 points. You can use an overlap than 1/2 of the FFT points to make smoother plots, but there is no increase of information for these higher overlaps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "nfft = 2048\n",
    "noverlap = nfft//2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another important parameter is the window function to use. We will use a Hanning window function named `signal.hann`. (This is actually the default window function used for the power spectrum, but we did not have to specify it in tht function call.)\n",
    "\n",
    "Enter the cell below to calcualate the spectrogram using the number of FFT points and overlap we specified and the Hanning window function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=[12, 6])\n",
    "pxx, freq, t, cax = ax.specgram(pipe_wave, Fs=pipe_info['fs'], NFFT=nfft)\n",
    "ax.set_xlabel('Time (s)')\n",
    "ax.set_ylabel('Frequency (Hz)')\n",
    "fig.colorbar(cax, label=\"log(PSD)\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can adjust some of the spectrogram parameters to show some of the sound characteristics more clearly. In the example below we use the variables `psdmin` and `psdmax` to set the minimum and maximum power spectral density (PSD) values for the color scale. We will also use `fmin` and `fmax` for the minimum and maximum frequencies displayed in the graph."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "psdmin = 10\n",
    "psdmax = 70\n",
    "fmin = 0\n",
    "fmax = 2000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=[12, 6])\n",
    "pxx, freq, t, cax = ax.specgram(pipe_wave, Fs=pipe_info['fs'], NFFT=nfft,\n",
    "                               vmin=psdmin, vmax=psdmax)\n",
    "ax.set_xlabel('Time (s)')\n",
    "ax.set_ylabel('Frequency (Hz)')\n",
    "ax.set_ylim([fmin, fmax])\n",
    "fig.colorbar(cax, label=\"log(PSD)\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that you can use a `ax.set_xlim` command (similar to the `set_ylim` command above) to display a smaller time interval."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Wave Plot and Spectrogram"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sometimes it is good to show both the waveform and the spectrogram of a sound in a combined graph with subplots."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(2, 1, sharex=True, figsize=[12, 10])\n",
    "ax[0].plot(pipe_times, pipe_wave, '-b')\n",
    "ax[0].set_ylabel('Relative Pressure')\n",
    "pxx, freq, t, cax = ax[1].specgram(pipe_wave, Fs=pipe_info['fs'], NFFT=nfft,\n",
    "                               vmin=psdmin, vmax=psdmax)\n",
    "ax[1].set_xlabel('Time (s)')\n",
    "ax[1].set_ylabel('Frequency (Hz)')\n",
    "ax[1].set_ylim([fmin, fmax])\n",
    "fig.colorbar(cax, label=\"log(PSD)\", orientation='horizontal', shrink=0.5)\n",
    "fig.tight_layout();"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SVPy",
   "language": "python",
   "name": "svpy"
  },
  "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.12.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}