connear_notebook_light.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "BKYqKPb26eeG"
   },
   "source": [
    "# *CoNNear*: A convolutional neural-network model of human cochlear mechanics and filter tuning for real-time applications\n",
    "\n",
    "This python notebook reproduces the evaluation results of the proposed CoNNear model. Please read through the prerequisites to install the correct packages."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "BKYqKPb26eeG"
   },
   "source": [
    "## Prerequisites\n",
    "\n",
    "- To run the trained CoNNear model, you should install numpy, scipy, keras and tensorflow, preferably through a Conda +3.6 environment. Open jupyter notebook to run this script."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 431
    },
    "colab_type": "code",
    "id": "WMcdw9aWCe0l",
    "outputId": "e6f8453e-4636-4f25-9a12-05e5aba07ce8",
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "import os\n",
    "#os.add_dll_directory(os.getcwd())\n",
    "\n",
    "import numpy as np\n",
    "from scipy import signal\n",
    "import scipy.signal as sp_sig\n",
    "import matplotlib.pyplot as plt\n",
    "import keras\n",
    "from keras.models import model_from_json\n",
    "from keras.utils import CustomObjectScope\n",
    "from keras.initializers import glorot_uniform\n",
    "\n",
    "from tlmodel.get_tl_vbm_and_oae import tl_vbm_and_oae\n",
    "from helper_ops import *\n",
    "\n",
    "json_file = open(\"connear/Gmodel.json\", \"r\")\n",
    "loaded_model_json = json_file.read()\n",
    "json_file.close()\n",
    "with CustomObjectScope({'GlorotUniform': glorot_uniform()}):    \n",
    "    connear = model_from_json(loaded_model_json)\n",
    "with CustomObjectScope({'GlorotUniform': glorot_uniform()}):\n",
    "     connear.load_weights(\"connear/Gmodel.h5\")\n",
    "connear.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "NqnSf4-8SUXE"
   },
   "source": [
    "Define parameters\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "PccC7PfWSNYB"
   },
   "outputs": [],
   "source": [
    "# Define model specific variables\n",
    "down_rate = 2\n",
    "fs = 20e3\n",
    "fs_tl = 100e3\n",
    "factor_fs = int(fs_tl / fs)\n",
    "p0 = 2e-5\n",
    "right_context = 256\n",
    "left_context = 256\n",
    "# load CFs \n",
    "CF = np.loadtxt('tlmodel/cf.txt')\n",
    "channels = CF.size "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "4otlWg_jS68e"
   },
   "source": [
    "## Click response\n",
    "Show the responses of the CoNNear model to a click stimulus.      \n",
    "**Notice that for all shown simulations the CoNNear model operates at 20kHz.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "yBGus8H8TFL5"
   },
   "outputs": [],
   "source": [
    "# Define the click stimulus\n",
    "dur = 128.0e-3 # for 2560 samples (2048 window length, 2x256 context)\n",
    "stim = np.zeros((1, int(dur * fs)))\n",
    "L = 70.0\n",
    "samples = dur * fs\n",
    "click_duration = 2 # 100 us click \n",
    "silence = 60 #samples in silence\n",
    "samples = int(samples - right_context - left_context)\n",
    "\n",
    "################ CoNNear ####################################\n",
    "stim = np.zeros((1, int(dur * fs)))\n",
    "stim[0, right_context + silence : right_context + silence + click_duration] = 2 * np.sqrt(2) * p0 * 10**(L/20)\n",
    "\n",
    "stim = np.expand_dims(stim, axis=2)\n",
    "connear_pred_click = connear.predict(stim.T, verbose=1)\n",
    "bmm_click_connear = connear_pred_click[0,:,:].T * 1e-6"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "5U7VOZU0R932"
   },
   "source": [
    "Plotting the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "D5_nN0QPRFtK"
   },
   "outputs": [],
   "source": [
    "################ Plots ######################################\n",
    "#Plot input stimulus\n",
    "plt.plot(stim[0,256:-256]), plt.xlim(0,2000)\n",
    "plt.show()\n",
    "\n",
    "# Plot the CoNNear response\n",
    "plt.imshow(bmm_click_connear, aspect='auto', cmap='jet')\n",
    "plt.xlim(0,2000), plt.clim(-4e-7,5e-7)\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cochlear Excitation Patterns\n",
    "Here, we plot the simulated RMS levels of basilar memberane (BM) displacement across CF for tone stimuli presented at SPLs between 0 and 90 dB SPL.\n",
    "\n",
    "**You can change the `f_tone` variable to simulate tone stimuli of different frequencies, say 500Hz, 1kHz, 2kHz, etc..**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f_tone = 1e3 # You can change this tone frequency to see how the excitation pattern changes with stimulus frequency\n",
    "dur = 102.4e-3 # for 2048 samples\n",
    "window_len = int(fs * dur)\n",
    "L = np.arange(0., 91.0, 10.) # SPLs from 0 to 90dB\n",
    "total_length = window_len + right_context + left_context #total length = 2560\n",
    "\n",
    "################ CoNNear ####################################\n",
    "t = np.arange(0., dur, 1./fs)\n",
    "hanlength = int(10e-3 * fs) # 10ms length hanning window\n",
    "stim_sin = np.sin(2 * np.pi * f_tone * t)\n",
    "han = signal.windows.hann(hanlength)\n",
    "stim_sin[:int(hanlength/2)] = stim_sin[:int(hanlength/2)] * han[:int(hanlength/2)]\n",
    "stim_sin[-int(hanlength/2):] = stim_sin[-int(hanlength/2):] * han[int(hanlength/2):]\n",
    "stim = np.zeros((len(L), total_length))\n",
    "for j in range(len(L)):\n",
    "    stim[j,right_context:window_len+right_context] = p0 * np.sqrt(2) * 10**(L[j]/20) * stim_sin\n",
    "    \n",
    "# prepare for feeding to the DNN\n",
    "stim = np.expand_dims(stim, axis=2)\n",
    "\n",
    "connear_pred_tone = connear.predict(stim, verbose=1)\n",
    "bmm_tone_connear = connear_pred_tone\n",
    "\n",
    "# Compute rms for each level\n",
    "cochlear_pred_tone_rms = np.vstack([rms(bmm_tone_connear[i]) for i in range(len(L))])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "################ Plots ######################################\n",
    "# Plot the RMS for CoNNear\n",
    "cftile=np.tile(CF, (len(L),1))\n",
    "plt.semilogx((cftile.T), 20.*np.log10(cochlear_pred_tone_rms.T))\n",
    "plt.xlim(0.25,8.), plt.grid(which='both'), \n",
    "plt.xticks(ticks=(0.25, 0.5, 1., 2., 4., 8.) , labels=(0.25, 0.5, 1., 2., 4., 8.))\n",
    "plt.ylim(-80, 20)\n",
    "plt.xlabel('CF (kHz)')\n",
    "plt.ylabel('RMS of y_bm (dB)')\n",
    "plt.title('CoNNear Predicted')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## QERB Plots\n",
    "\n",
    "Next, the level dependent tuning properties of the cochlear filters (QERB) are shown. We chose click stimuli at three levels 0dB, 40dB and 70dB peSPL, computed the CoNNEar cochlear BM response and computed the QERB."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "L = [0., 40., 70.] # We will plot it for three SPLs\n",
    "    \n",
    "################ CoNNear ####################################\n",
    "stim_con = np.zeros(((len(L),int(dur * fs)+right_context+left_context,1)))\n",
    "QERB_connear = np.zeros(((len(L),channels)))\n",
    "\n",
    "for i in range (len(L)):\n",
    "    stim_con[i, right_context + silence : right_context + silence + click_duration, 0] = 2 * np.sqrt(2) * p0 * 10**(L[i]/20)\n",
    "    \n",
    "#Get CoNNear outputs\n",
    "con_predicted = connear.predict(stim_con)\n",
    "\n",
    "for i in range (len(L)):\n",
    "    QERB_connear[i,:] = QERB_calculation(con_predicted[i, :, :].T, CF*1e3, fs)\n",
    "\n",
    "################ Plots ######################################\n",
    "# Plot QERB of CoNNear model\n",
    "plt.semilogx(CF[0::5], (QERB_connear[0,0::5]),':gs', label='0dB')\n",
    "plt.semilogx(CF[0::5], (QERB_connear[1,0::5]),'r', label='40dB')\n",
    "plt.semilogx(CF[0::5], (QERB_connear[2,0::5]),':rs', label='70dB')\n",
    "plt.xlim(0.25,8.), plt.grid(which='both'), \n",
    "plt.xticks(ticks=(0.25, 0.5, 1., 2., 4., 8.) , labels=(0.25, 0.5, 1., 2., 4., 8.))\n",
    "plt.yticks(ticks=(5, 10, 15, 20) , labels=(5, 10, 15, 20))\n",
    "plt.xlabel('CF (kHz)')\n",
    "plt.ylim(2,20)\n",
    "plt.ylabel('QERB')\n",
    "plt.title('CoNNear Predicted')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Speech Input\n",
    "Here, a sentence from the Dutch speech matrix (unseen during training) will be input to both the TL and the CoNNear models. By adapting fragment_length parameters, various input lengths can be compared and visualised. \n",
    "\n",
    "**Notice that this part is computationally more expensive with a higher fragment_duration, both in terms of memory and time.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#load in speechfile\n",
    "wavfile = 'dutch_sentence.wav'\n",
    "signal_wav, fs_signal = wavfile_read(wavfile)\n",
    "signalr = sp_sig.resample_poly(signal_wav, fs_tl, fs_signal)\n",
    "\n",
    "L = np.array([70]) #sound-level of 70 dB SPL\n",
    "stim_full = np.zeros((len(L), signalr.size))\n",
    "\n",
    "for j in range(len(L)):\n",
    "    stim_full[j, :] = p0 * 10**(L[j]/20) * signalr/rms(signalr)\n",
    "    \n",
    "fragment_length = 12345 #define fragment length (max 40000 for the included wav-file)\n",
    "stim_length_init = factor_fs*(fragment_length+right_context+left_context)\n",
    "stim_length = stim_length_init\n",
    "\n",
    "#adapt fragment duration if no multiple of 16 (due to the CNN character of CoNNear, we need multiples of 16)\n",
    "zero_pad = fragment_length%16\n",
    "zeros = 0\n",
    "if zero_pad != 0: \n",
    "    zeros = 16-zero_pad\n",
    "    stim_length = factor_fs*(fragment_length+right_context+left_context+zeros)\n",
    "        \n",
    "################ CoNNear ####################################\n",
    "stim = np.zeros((len(L), int(stim_length)))\n",
    "stimrange = range(0, stim_length_init)\n",
    "stim[:,stimrange] = stim_full[:,0:stim_length_init]\n",
    "stim=sp_sig.resample_poly(stim, fs, fs_tl, axis=1)\n",
    "stim=np.expand_dims(stim, axis=2)\n",
    "\n",
    "tl_pred = connear.predict(stim)\n",
    "tl_pred = tl_pred[0, :, :].T\n",
    "\n",
    "################ Plots ######################################\n",
    "fig, axarr = plt.subplots(2, sharex=True)\n",
    "axarr[0].set_ylim(-0.35, 0.35)\n",
    "axarr[0].plot(stim[0,(left_context):-(right_context),0])\n",
    "axarr[0].set_title('Segment of Audio Input')\n",
    "\n",
    "cax2 = axarr[1].imshow(tl_pred, cmap='bwr',aspect='auto', vmin=-0.5, vmax=0.5)  \n",
    "axarr[1].set_title('Output of CoNNear')\n",
    "axarr[1].set(ylabel='Center Frequency')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## DPOAE Plots\n",
    "The frequency response of the 12-kHz CF channel is evaluated as a proxy for the otoacoustic emissions recorded in the ear-canal. Frequency responses of model simulations are shown in response to two pure tones of $f_{1,2}$ of 2.0 and 2.4 kHz.\n",
    "The most pronounced distortion product in humans occurs at $2f_1 - f_2$ (1.6 kHz).\n",
    "\n",
    "**Notice that this part is computationally expensive (up to 20 mins), both in terms of memory and time.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create the stimulus\n",
    "L = [70.0]\n",
    "f1 = 2000.\n",
    "f2 = 1.2 * f1\n",
    "L2 = 50.\n",
    "L1 = 39 + 0.4 * L2 # scissors paradigm\n",
    "print (\"The tone frequencies are \" + str(f1) + \" and \" + str(f2))\n",
    "print (\"with levels \" + str(L1) + \" and \" + str(L2))\n",
    "trailing_silence = 0.\n",
    "\n",
    "# Here we will a pick stimulus longer than 2048 samples to get a better FFT\n",
    "# We will prepare it with fs_tl sampling frequency\n",
    "dur_sin_samples = np.lcm(int(f1), int(f2))\n",
    "min_duration = 0.25 # in seconds\n",
    "if dur_sin_samples < (min_duration * fs_tl):\n",
    "    dur_sin = (((min_duration * fs_tl) // dur_sin_samples) + 1) * dur_sin_samples\n",
    "else:\n",
    "    dur_sin = dur_sin_samples\n",
    "dur_sin = (dur_sin / fs_tl) \n",
    "\n",
    "t = np.arange(0., dur_sin, 1./fs_tl)\n",
    "\n",
    "hanlength = int(10e-3 * fs_tl) # 10ms length hanning window\n",
    "#f1\n",
    "stim_sin1 = np.sin(2 * np.pi * f1 * t)\n",
    "han = signal.windows.hann(hanlength)\n",
    "stim_sin1 = p0 * np.sqrt(2) * 10**(L1/20) * stim_sin1\n",
    "#f2\n",
    "stim_sin2 = np.sin(2 * np.pi * f2 * t)\n",
    "stim_sin2 = p0 * np.sqrt(2) * 10**(L2/20) * stim_sin2\n",
    "\n",
    "stim_sin = stim_sin1 + stim_sin2\n",
    "\n",
    "total_length =  int(trailing_silence * fs_tl) + len(stim_sin)\n",
    "stim = np.zeros((1, int(total_length)))\n",
    "stimrange = range(int(trailing_silence * fs_tl), int(trailing_silence * fs_tl) + len(stim_sin))\n",
    "stim[0, stimrange] = stim_sin\n",
    "\n",
    "################# CoNNear ####################################\n",
    "# prepare for feeding to the CoNNear model\n",
    "# first resample it to fs\n",
    "shift_stim = 1\n",
    "stim = stim[:, :]\n",
    "stim = signal.decimate(stim, factor_fs, axis=1)\n",
    "stim_1 = np.array(stim[0,:])\n",
    "# window the signal into chunks of 2560 samples to be fed to the CoNNer model\n",
    "stim = slice_1dsignal(stim_1, 2048, shift_stim, 256, left_context=256, right_context=256) \n",
    "connear_out_chunks = connear.predict(stim,verbose=1)\n",
    "# undo the windowing to get back the full response\n",
    "connear_out_full = undo_window (connear_out_chunks, 2048, shift_stim, ignore_first_set=0)\n",
    "connear_out_full = connear_out_full[:,:stim_1.shape[0],:] * 1e-6\n",
    "\n",
    "##############################################################\n",
    "f_cf = 12000.\n",
    "tone_index, tone_cf = min(enumerate(CF*1000), key=lambda x: abs( x [1]- f_cf))\n",
    "print(\"CF nearest to \" + str(f_cf) + \" is \" + str(CF[tone_index]*1000))\n",
    "scale_val = (p0* np.sqrt(2))\n",
    "\n",
    "################ Plots ######################################\n",
    "# Plot the DPOAE CoNNear\n",
    "connear_dpoae, nfft_connear = get_dpoae(connear_out_full, cf_location=tone_index)\n",
    "    \n",
    "freq_bins_connear = np.linspace(0, fs, num = nfft_connear)\n",
    "plt.semilogx(freq_bins_connear[:int(nfft_connear/2)]/1000, 20 * np.log10(connear_dpoae/scale_val))\n",
    "plt.title(\"DPOAEs - CoNNear\")\n",
    "plt.xlabel('Frequency [kHz]'), plt.ylabel('Magnitudes dbSPL'), plt.xlim((0.25, 8)), plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## RMS error between the excitation patterns\n",
    "\n",
    "This is not covered in the light version of the notebook."
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "name": "connear_notebook.ipynb",
   "provenance": [],
   "toc_visible": true
  },
  "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.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}