gin commit from DESKTOP-CL7Q894

New files: 13

ERBFilterBank.m

@@ -0,0 +1,60 @@
+function output = ERBFilterBank(x, fcoefs)
+% function output = ERBFilterBank(x, fcoefs)
+% Process an input waveform with a gammatone filter bank. This function 
+% takes a single sound vector, and returns an array of filter outputs, one 
+% channel per row.
+%
+% The fcoefs parameter, which completely specifies the Gammatone filterbank,
+% should be designed with the MakeERBFilters function.  If it is omitted,
+% the filter coefficients are computed for you assuming a 22050Hz sampling
+% rate and 64 filters regularly spaced on an ERB scale from fs/2 down to 100Hz.
+%
+
+% Malcolm Slaney @ Interval, June 11, 1998.
+% (c) 1998 Interval Research Corporation  
+% Thanks to Alain de Cheveigne' for his suggestions and improvements.
+
+if nargin < 1
+	error('Syntax: output_array = ERBFilterBank(input_vector[, fcoefs]);');
+end
+
+if nargin < 2
+	fcoefs = MakeERBFilters(22050,64,100);
+end
+
+if size(fcoefs,2) ~= 10
+	error('fcoefs parameter passed to ERBFilterBank is the wrong size.');
+end
+
+if size(x,2) < size(x,1)
+	x = x';
+end
+
+A0  = fcoefs(:,1);
+A11 = fcoefs(:,2);
+A12 = fcoefs(:,3);
+A13 = fcoefs(:,4);
+A14 = fcoefs(:,5);
+A2  = fcoefs(:,6);
+B0  = fcoefs(:,7);
+B1  = fcoefs(:,8);
+B2  = fcoefs(:,9);
+gain= fcoefs(:,10);	
+
+output = zeros(size(gain,1), length(x));
+for chan = 1: size(gain,1)
+	y1=filter([A0(chan)/gain(chan) A11(chan)/gain(chan) ...
+		   A2(chan)/gain(chan)], ...
+				[B0(chan) B1(chan) B2(chan)], x);
+	y2=filter([A0(chan) A12(chan) A2(chan)], ...
+				[B0(chan) B1(chan) B2(chan)], y1);
+	y3=filter([A0(chan) A13(chan) A2(chan)], ...
+				[B0(chan) B1(chan) B2(chan)], y2);
+	y4=filter([A0(chan) A14(chan) A2(chan)], ...
+				[B0(chan) B1(chan) B2(chan)], y3);
+	output(chan, :) = y4;
+end
+
+if 0
+	semilogx((0:(length(x)-1))*(fs/length(x)),20*log10(abs(fft(output))));
+end

Get_Spike_Counts.ipynb

@@ -0,0 +1,464 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Abi"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "abis = np.asarray([trial['abi'] for trial in abifile.trials])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "spiketrains = np.empty((len(abifile.events), ), dtype='object')\n",
+    "for kt, tevents in enumerate(abifile.events):\n",
+    "    spiketimes = np.asarray([t for t, unit in tevents])\n",
+    "    spiketrains[kt] = spiketimes / 10000\n",
+    "order = np.argsort(abis)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "scrolled": true
+   },
+   "outputs": [],
+   "source": [
+    "abispikecounts = np.zeros((abifile.params['abi']['range'].size, abifile.params['reps']))\n",
+    "for kabi, abi in enumerate(abifile.params['abi']['range']):\n",
+    "    for krep, spiketimes in enumerate(spiketrains[abis == abi]):\n",
+    "        abispikecounts[kabi, krep] = within(spiketimes, bffile.params['delay'], bffile.params['delay'] + \\\n",
+    "                                        bffile.params['dur']*2)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## ILD"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ilds = np.asarray([trial['iid'] for trial in ildfile.trials])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "spiketrains = np.empty((len(ildfile.events), ), dtype='object')\n",
+    "for kt, tevents in enumerate(ildfile.events):\n",
+    "    spiketimes = np.asarray([t for t, unit in tevents])\n",
+    "    spiketrains[kt] = spiketimes / 10000\n",
+    "order = np.argsort(ilds)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "scrolled": true
+   },
+   "outputs": [],
+   "source": [
+    "ildspikecounts = np.zeros((ildfile.params['iid']['range'].size, ildfile.params['reps']))\n",
+    "for kabi, iid in enumerate(ildfile.params['iid']['range']):\n",
+    "    for krep, spiketimes in enumerate(spiketrains[ilds == iid]):\n",
+    "        ildspikecounts[kabi, krep] = within(spiketimes, bffile.params['delay'], bffile.params['delay'] + \\\n",
+    "                                        bffile.params['dur']*2)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ildx = ildfile.params['iid']['range']\n",
+    "ildx = np.round(ildx, decimals=2)\n",
+    "ildy = ildspikecounts.mean(axis=1)\n",
+    "ildy = np.round(ildy, decimals=2)\n",
+    "ildz = sps.stats.sem(ildspikecounts, axis=1)\n",
+    "ildz[ildz==0] = 0.01\n",
+    "so = [np.max(ildy), 1, 0, 1, np.min(ildy)]\n",
+    "go = [np.min(ildy), np.max(ildy), ildx[np.argmax(ildy)], 10]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "try:\n",
+    "    gopt, gcov = spo.curve_fit(Gaussian, ildx, ildy, go, sigma = ildz, absolute_sigma = True)\n",
+    "    gcorr = np.corrcoef(ildy, Gaussian(ildx, *gopt))\n",
+    "    Rsqg = gcorr[0,1]**2\n",
+    "except RuntimeError:\n",
+    "    print('No good Gaussian fit')\n",
+    "    Rsqg = np.nan\n",
+    "try:\n",
+    "    sopt, scov = spo.curve_fit(Sigmoidal, ildx, ildy, so, sigma = ildz, absolute_sigma = True)\n",
+    "    scorr = np.corrcoef(ildy, Sigmoidal(ildx, *sopt))\n",
+    "    Rsqs = scorr[0,1]**2\n",
+    "except RuntimeError:\n",
+    "    print('no good Sigmoidal fit')\n",
+    "    Rsqs = np.nan"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "peak = np.where(ildy == ildy.max())\n",
+    "ildhalfheight = np.where(ildy == find_nearest(ildy,(ildy.max() - ildy.min())/2 ) )\n",
+    "if ildhalfheight[0].size == 0:\n",
+    "    ildhalfwidth = np.empty([0])\n",
+    "    ildhalfwidth = np.nan    \n",
+    "elif ildhalfheight[0].size == 1:\n",
+    "    ildhalfwidth = np.abs(ildx[peak] - ildx[ildhalfheight])\n",
+    "elif ildhalfheight[0].size == 2: \n",
+    "    ildwidth = ildhalfheight[0][1] - ildhalfheight[0][0]\n",
+    "best_ild = ildx[peak]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "##### BF"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "bfs = np.asarray([trial['stim'] for trial in bffile.trials])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "bfspiketrains = np.empty((len(bffile.events), ), dtype='object')\n",
+    "for kt, tevents in enumerate(bffile.events):\n",
+    "    spiketimes = np.asarray([t for t, unit in tevents])\n",
+    "    bfspiketrains[kt] = spiketimes / 10000\n",
+    "bforder = np.argsort(bfs)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "scrolled": true
+   },
+   "outputs": [],
+   "source": [
+    "bfspikecounts = np.zeros((bffile.params['bf']['range'].size, bffile.params['reps']))\n",
+    "\n",
+    "for kabi, stim in enumerate(bffile.params['bf']['range']):\n",
+    "    for krep, spiketimes in enumerate(bfspiketrains[bfs == stim]):\n",
+    "        bfspikecounts[kabi, krep] = within(spiketimes, bffile.params['delay'], bffile.params['delay'] + \\\n",
+    "                                        bffile.params['dur']*2)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "bfx = bffile.params['bf']['range']\n",
+    "bfy = bfspikecounts.mean(axis=1)\n",
+    "bfse = sps.stats.sem(bfspikecounts, axis=1)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "peak_idx = np.argmax(bfy)\n",
+    "halfmax_y = (np.max(bfy) - np.min(bfy))/2 + np.min(bfy)\n",
+    "\n",
+    "lup = np.where(bfy[:peak_idx] > halfmax_y)[0]\n",
+    "ldown = np.where(bfy[:peak_idx] < halfmax_y)[0]\n",
+    "rup = np.where(bfy[peak_idx+1:] > halfmax_y)[0]\n",
+    "rdown = np.where(bfy[peak_idx+1:] < halfmax_y)[0]\n",
+    "pm = np.array([-1,1])\n",
+    "\n",
+    "\n",
+    "if lup.size > 0 and ldown.size > 0:\n",
+    "    if between(halfmax_y, bfy[lup[0]] +  pm*bfse[lup[0]]):\n",
+    "        l = lup[0]\n",
+    "    elif between(halfmax_y, bfy[ldown[0]] +  pm*bfse[ldown[0]]):\n",
+    "        l = ldown[-1] \n",
+    "    elif lup[0] == ldown[-1] + 1:\n",
+    "        l = lup[0] \n",
+    "    else:\n",
+    "        l = lup[0]\n",
+    "elif lup.size > 0:\n",
+    "    l = lup[0]\n",
+    "elif ldown.size > 0: \n",
+    "    l = ldown[-1] \n",
+    "else:\n",
+    "    l = peak_idx\n",
+    "\n",
+    "if rup.size > 0 and rdown.size > 0:\n",
+    "    if between(halfmax_y, bfy[rup[-1] + peak_idx + 1] +  pm*bfse[rup[-1] + peak_idx + 1]):\n",
+    "        r = rup[-1] + peak_idx + 1\n",
+    "    elif between(halfmax_y, bfy[rdown[0] + peak_idx + 1] +  pm*bfse[rdown[0] + peak_idx + 1]):\n",
+    "        r = rdown[0] + peak_idx +1\n",
+    "    elif rup[-1] == rdown[0] - 1:\n",
+    "        r = rup[-1] + peak_idx + 1\n",
+    "    else:\n",
+    "        r = rup[-1] + peak_idx + 1\n",
+    "elif rup.size > 0:\n",
+    "    r = rup[-1] + peak_idx + 1\n",
+    "elif rdown.size > 0: \n",
+    "    r = rdown[0] + peak_idx +1\n",
+    "else:\n",
+    "    r = peak_idx   \n",
+    "\n",
+    "    \n",
+    "freq_lo, lofr = bfx[l], bfy[l]\n",
+    "freq_hi, hifr = bfx[r], bfy[r]\n",
+    "best_freq = np.mean([freq_lo, freq_hi])"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "**ITD**"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "itds = np.asarray([trial['itd'] for trial in itdfile.trials])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "itdspiketrains = np.empty((len(itdfile.events), ), dtype='object')\n",
+    "for kt, tevents in enumerate(itdfile.events):\n",
+    "    spiketimes = np.asarray([t for t, unit in tevents])\n",
+    "    itdspiketrains[kt] = spiketimes / 10000\n",
+    "itdorder = np.argsort(itds)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "scrolled": true
+   },
+   "outputs": [],
+   "source": [
+    "itdspikecounts = np.zeros((itdfile.params['itd']['range'].size, itdfile.params['reps']))\n",
+    "for kabi, stim in enumerate(itdfile.params['itd']['range']):\n",
+    "    for krep, spiketimes in enumerate(itdspiketrains[itds == stim]):\n",
+    "        itdspikecounts[kabi, krep] = within(spiketimes, itdfile.params['delay'], itdfile.params['delay'] + \\\n",
+    "                                        itdfile.params['dur']*2)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "itdx = itdfile.params['itd']['range']\n",
+    "itdy = itdspikecounts.mean(axis=1)\n",
+    "itdz = sps.stats.sem(itdspikecounts, axis=1)\n",
+    "itdz[itdz==0] = 0.01\n",
+    "H = lambda n: (bfx[n+1]-bfx[n]) * ((halfmax_y-bfy[n]) / (bfy[n+1]-bfy[n])) + bfx[n]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "scrolled": true
+   },
+   "outputs": [],
+   "source": [
+    "peak_idx = np.argmax(itdy)\n",
+    "halfmax_y = (np.max(itdy) - np.min(itdy))/2 + np.min(itdy)\n",
+    "#halfmax_y = np.mean(itdy) + np.std(itdy)\n",
+    "if np.where(itdy[:peak_idx] >= halfmax_y)[0].size > 0:\n",
+    "    l = np.where(itdy[:peak_idx] >= halfmax_y)[0][0]\n",
+    "else: #if no values besides peak is greater than halfmax_y\n",
+    "    l = peak_idx - 1          \n",
+    "if np.where(itdy[peak_idx+1:] >= halfmax_y)[0].size > 0: #peak__idx + 1 due to indexing issues\n",
+    "    r = np.where(itdy[peak_idx+1:] >= halfmax_y)[0][-1] + peak_idx +1\n",
+    "    if r == itdx.argmax():\n",
+    "        r = r-1\n",
+    "else:\n",
+    "    r = peak_idx + 1 "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "scrolled": true
+   },
+   "outputs": [],
+   "source": [
+    "if np.size(spg.find_peaks(itdy, height = np.max(itdy)*.8, distance = 150/10, prominence = 1)[0]) == 3:\n",
+    "    [peak, locs] = spg.find_peaks(itdy, height = np.max(itdy)*.8, distance = 150/10, prominence = 1)\n",
+    "    mainpk = np.arange(locs['left_bases'][1],locs['right_bases'][1])\n",
+    "    po = np.array([np.min(itdy), locs['peak_heights'][1]- np.min(itdy[mainpk]), itdx[peak[1]], 50],dtype = object)\n",
+    "elif np.size(spg.find_peaks(itdy, height = np.max(itdy)*.8, distance = 150/10, prominence = 1.4)[0]) == 1: \n",
+    "    [peak, locs] = spg.find_peaks(itdy, height = np.max(itdy)*.8, distance = 150/10, prominence = 1.4)\n",
+    "    mainpk = np.arange(locs['left_bases'][0],locs['right_bases'][0])\n",
+    "    po = np.array([np.min(itdy), locs['peak_heights'][0]- np.min(itdy[mainpk]), itdx[peak[0]], 50],dtype = object)\n",
+    "elif np.size(spg.find_peaks(itdy, height = np.max(itdy)*.8, distance = 150/10, prominence = 1.4)[0]) > 0:\n",
+    "    [peak, locs] = spg.find_peaks(itdy, height = np.max(itdy)*.8, distance = 150/10, prominence = 1.4)\n",
+    "    ppk = locs['peak_heights'].argmax()\n",
+    "    mainpk = np.arange(locs['left_bases'][ppk],locs['right_bases'][ppk])\n",
+    "    po = np.array([np.min(itdy), locs['peak_heights'][ppk]- np.min(itdy[mainpk]), itdx[peak[ppk]], 50],dtype = object)\n",
+    "if np.size(spg.find_peaks(itdy, height = np.max(itdy)*.8, distance = 150/10, prominence = 1.4)[0]) > 0:\n",
+    "    try:\n",
+    "        bnds = [np.min(itdy) - itdz[np.argmin(itdy)], np.max(itdy[mainpk]) - np.min(itdy[mainpk]) - itdz[np.argmin(itdy[mainpk])], -550, 1], \\\n",
+    "               [np.min(itdy) + itdz[np.argmin(itdy)], np.max(itdy[mainpk]) - np.min(itdy[mainpk]) + itdz[np.argmax(itdy[mainpk])], 550, 300]\n",
+    "        ITDopt, ITDcov = spo.curve_fit(Gaussian, itdx[mainpk], itdy[mainpk], po, sigma = itdz[mainpk],\n",
+    "                                       bounds = bnds, absolute_sigma = True)\n",
+    "        ITDcorr = np.corrcoef(itdy[mainpk], Gaussian(itdx[mainpk], *ITDopt))\n",
+    "        ITDRsq = ITDcorr[0,1]**2\n",
+    "    except RuntimeError:\n",
+    "        print('no good Gaussian fit')\n",
+    "        ITDRsq = np.nan\n",
+    "else: \n",
+    "    ITDRsq = np.nan"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "scrolled": true
+   },
+   "outputs": [],
+   "source": [
+    "H = lambda n: (itdx[n+1]-itdx[n]) * ((halfmax_y-itdy[n]) / (itdy[n+1]-itdy[n])) + itdx[n]\n",
+    "if ITDRsq >.5:\n",
+    "    gITD = Gaussian(itdx, *ITDopt)\n",
+    "    xint = np.arange(itdx[0],itdx[-1])\n",
+    "    gint = np.interp(np.arange(itdx[0],itdx[-1]), itdx, gITD)\n",
+    "    gpeak_idx = np.argmax(gint)\n",
+    "    ghalfpeak = (np.max(gint) - np.min(gint))/2 + np.min(gint)\n",
+    "    best_itd = xint[gpeak_idx]\n",
+    "    ghalfwidth = np.abs(xint[np.abs(gint-ghalfpeak).argmin()] - xint[gpeak_idx])\n",
+    "    ITD_lo = xint[gpeak_idx] - ghalfwidth\n",
+    "    ITD_hi = xint[gpeak_idx] + ghalfwidth\n",
+    "    #display([np.mean((H(l), H(r))).round(2),xint[gpeak_idx]], \\\n",
+    "            #[H(l).round(2), xint[gpeak_idx] - ghalfwidth], [H(r).round(2),xint[gpeak_idx] + ghalfwidth])\n",
+    "else:\n",
+    "    best_itd = np.mean((H(l), H(r)))\n",
+    "    ITD_lo = H(l)\n",
+    "    ITD_hi = H(r)\n",
+    "    #display([np.mean((H(l), H(r)))], [H(l)], [H(r)])\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "scrolled": true
+   },
+   "outputs": [],
+   "source": [
+    "if ~np.isnan(ITDRsq):\n",
+    "    print('best ITD:')\n",
+    "    display([best_itd.round(), ITD_lo.round(), ITD_hi.round(), ITDRsq.round(2)])\n",
+    "    print('best freq:')\n",
+    "    display([best_freq.round(),freq_lo.round(),freq_hi.round()])\n",
+    "    print('best ILD:')\n",
+    "    display(best_ild)\n",
+    "    print('Gaussian or Sigmoidal')\n",
+    "    display(np.round([Rsqg, Rsqs], decimals = 2))\n",
+    "else:\n",
+    "    print('best ITD:')\n",
+    "    display([best_itd.round(), ITD_lo.round(), ITD_hi.round(), ITDRsq])\n",
+    "    print('best freq:')\n",
+    "    display([best_freq.round(),freq_lo.round(),freq_hi.round()])\n",
+    "    print('best ILD:')\n",
+    "    display(best_ild)\n",
+    "    print('Gaussian or Sigmoidal')\n",
+    "    display(np.round([Rsqg, Rsqs], decimals = 2))"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "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.9.12"
+  },
+  "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
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}

Read_XDPhys_Raw_Data.ipynb

@@ -0,0 +1,438 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import pyxdphys\n",
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "import scipy.stats as sps\n",
+    "import scipy.optimize as spo\n",
+    "import scipy.signal as spg\n",
+    "import matplotlib as mpl\n",
+    "import matplotlib.pyplot as plt"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# matplotlib settings\n",
+    "mpl.rcParams['mathtext.default'] = 'regular'\n",
+    "\n",
+    "# Reproducible SVG files (including xml ids):\n",
+    "mpl.rcParams['svg.hashsalt'] = '123'\n",
+    "mpl.rcParams['savefig.dpi'] = 300\n",
+    "\n",
+    "# Font sizes:\n",
+    "mpl.rcParams['axes.labelsize'] =  18 # default: 'medium' == 10\n",
+    "mpl.rcParams['xtick.labelsize'] = 18 # default: 'medium' == 10\n",
+    "mpl.rcParams['ytick.labelsize'] = 18 # default: 'medium' == 10\n",
+    "mpl.rcParams['legend.fontsize'] = 18\n",
+    "\n",
+    "# Other texts:\n",
+    "in_panel_fontsize = 18\n",
+    "bc_indicator_size = in_panel_fontsize"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def find_nearest(array, value):\n",
+    "    array = np.asarray(array)\n",
+    "    idx = (np.abs(array - value)).argmin()\n",
+    "    return [idx]\n",
+    "\n",
+    "def Gaussian(t, ga, gb, gc, gd):\n",
+    "    return ga + np.abs(gb) * np.exp( - ( (t-gc)**2 / gd**2) )\n",
+    "\n",
+    "def Sigmoidal(t, sa, sb, sc, sd, se):\n",
+    "    return sa / (1 + np.exp(-sb*(t - sc)/sd)) + se\n",
+    "\n",
+    "def within(tval, tmin, tmax):\n",
+    "    return np.sum((tval > tmin) & (tval < tmax))\n",
+    "\n",
+    "def between(val, array):\n",
+    "    return array[0] <= val <= array[1]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {
+    "scrolled": false
+   },
+   "outputs": [],
+   "source": [
+    "## Select neuron. d = date, m = owl.neuron\n",
+    "d = '2021-0929'\n",
+    "m =  '44.02'\n",
+    "n = '0' + m\n",
+    "\n",
+    "#alt, in cases where abi was recorded using using monaural and binuaral (L,R,B):\n",
+    "#abifile = pyxdphys.XDdata(fr'F:\\xdphys\\{d}\\{m}\\{n}.0.1.abi')\n",
+    "abifile  = pyxdphys.XDdata(fr'F:\\xdphys\\{d}\\{m}\\{n}.0.abi')\n",
+    "ildfile  = pyxdphys.XDdata(fr'F:\\xdphys\\{d}\\{m}\\{n}.1.iid')\n",
+    "bcfile   = pyxdphys.XDdata(fr'F:\\xdphys\\{d}\\{m}\\{n}.2.gen')\n",
+    "itdfile  = pyxdphys.XDdata(fr'F:\\xdphys\\{d}\\{m}\\{n}.3.itd')\n",
+    "bffile   = pyxdphys.XDdata(fr'F:\\xdphys\\{d}\\{m}\\{n}.4.bf')\n",
+    "#sptfile  = pyxdphys.XDdata(r'D:\\xdphys\\2021-0811\\046.00\\046.00.5.iid')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {
+    "scrolled": true
+   },
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "best ITD:\n"
+     ]
+    },
+    {
+     "data": {
+      "text/plain": [
+       "[20.0, -21.0, 61.0, 0.99]"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    },
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "best freq:\n"
+     ]
+    },
+    {
+     "data": {
+      "text/plain": [
+       "[4900.0, 4000.0, 5800.0]"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    },
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "best ILD:\n"
+     ]
+    },
+    {
+     "data": {
+      "text/plain": [
+       "array([4.])"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    },
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Gaussian or Sigmoidal\n"
+     ]
+    },
+    {
+     "data": {
+      "text/plain": [
+       "array([0.99, 0.08])"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "##Running this will do all analysis automatically. The rest of this script is in case there are issues,\n",
+    "## or to see rasterplots/plots\n",
+    "%run ./Get_Spike_Counts.ipynb"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#run next cell if best ITD isn't found. This usually fixes the issue. xlsx notes when this is used"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "siteITD = 0\n",
+    "[peak, locs] = spg.find_peaks(itdy, height = np.max(itdy)*.8, distance = 100/10, prominence = 1.)\n",
+    "\n",
+    "tt = find_nearest(itdx[peak], siteITD)\n",
+    "\n",
+    "mainpk = np.arange(locs['left_bases'][tt],locs['right_bases'][tt])\n",
+    "po = np.array([np.min(itdy), locs['peak_heights'][tt]- np.min(itdy[mainpk]), itdx[peak[tt]], 50],dtype=object)\n",
+    "\n",
+    "try:\n",
+    "    bnds = [[np.min(itdy) - itdz[np.argmin(itdy)], np.max(itdy[mainpk]) - np.min(itdy[mainpk]) - itdz[np.argmin(itdy[mainpk])], itdx[mainpk].min() - 100, 1], \\\n",
+    "           [np.min(itdy) + itdz[np.argmin(itdy)], np.max(itdy[mainpk]) - np.min(itdy[mainpk]) + itdz[np.argmax(itdy[mainpk])], itdx[mainpk].max() + 100, 300]]\n",
+    "    ITDopt, ITDcov = spo.curve_fit(Gaussian, itdx[mainpk], itdy[mainpk], po, sigma = itdz[mainpk],\n",
+    "                                   bounds = bnds, absolute_sigma = True)\n",
+    "    ITDcorr = np.corrcoef(itdy[mainpk], Gaussian(itdx[mainpk], *ITDopt))\n",
+    "    ITDRsq = ITDcorr[0,1]**2\n",
+    "except RuntimeError:\n",
+    "    print('no good Gaussian fit')\n",
+    "    ITDRsq = np.nan\n",
+    "    \n",
+    "H = lambda n: (itdx[n+1]-itdx[n]) * ((halfmax_y-itdy[n]) / (itdy[n+1]-itdy[n])) + itdx[n]\n",
+    "if ITDRsq >.25:\n",
+    "    gITD = Gaussian(itdx, *ITDopt)\n",
+    "    xint = np.arange(itdx[0],itdx[-1])\n",
+    "    gint = np.interp(np.arange(itdx[0],itdx[-1]), itdx, gITD)\n",
+    "    gpeak_idx = np.argmax(gint)\n",
+    "    ghalfpeak = (np.max(gint) - np.min(gint))/2 + np.min(gint)\n",
+    "    best_itd = xint[gpeak_idx]\n",
+    "    ghalfwidth = np.abs(xint[np.abs(gint-ghalfpeak).argmin()] - xint[gpeak_idx])\n",
+    "    ITD_lo = xint[gpeak_idx] - ghalfwidth\n",
+    "    ITD_hi = xint[gpeak_idx] + ghalfwidth\n",
+    "    display([np.mean((H(l), H(r))).round(2),xint[gpeak_idx]], \\\n",
+    "            [H(l).round(2), xint[gpeak_idx] - ghalfwidth], [H(r).round(2),xint[gpeak_idx] + ghalfwidth])\n",
+    "else:\n",
+    "    best_itd = np.mean((H(l), H(r)))\n",
+    "    ITD_lo = H(l)\n",
+    "    ITD_hi = H(r)\n",
+    "    display([np.mean((H(l), H(r)))], [H(l)], [H(r)])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#Best ITD Graph\n",
+    "fig, ax = plt.subplots(figsize = [5,4])\n",
+    "ax.errorbar(itdfile.params['itd']['range'], itdspikecounts.mean(axis=1), sps.stats.sem(itdspikecounts, axis=1))\n",
+    "ax.plot(itdx, Gaussian(itdx, *ITDopt))\n",
+    "plt.ylabel('Spike Count')\n",
+    "plt.xlabel('ITD (us)')\n",
+    "\n",
+    "#fig.savefig('ExampleITD.png')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# ITD raster plot\n",
+    "fig, ax = plt.subplots()\n",
+    "ax.eventplot(itdspiketrains[itdorder], color='k')\n",
+    "ax.set_yticks(np.arange(itdfile.params['itd']['range'].size)[::4] \\\n",
+    "              * itdfile.params['reps'] + itdfile.params['reps']/2)\n",
+    "ax.set_yticklabels(itdfile.params['itd']['range'][::4])\n",
+    "\n",
+    "ax.axvspan(itdfile.params['delay'], itdfile.params['delay']+itdfile.params['dur']+10, color='yellow')\n",
+    "ax.axvspan(itdfile.params['delay'], itdfile.params['delay']+ 15, color='pink')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#run the next four cells to change bffile analysis window to just the sound duration\n",
+    "bfs = np.asarray([trial['stim'] for trial in bffile.trials])\n",
+    "bfspiketrains = np.empty((len(bffile.events), ), dtype='object')\n",
+    "for kt, tevents in enumerate(bffile.events):\n",
+    "    spiketimes = np.asarray([t for t, unit in tevents])\n",
+    "    bfspiketrains[kt] = spiketimes / 10000\n",
+    "bforder = np.argsort(bfs)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "bfspikecounts = np.zeros((bffile.params['bf']['range'].size, bffile.params['reps']))\n",
+    "\n",
+    "for kabi, stim in enumerate(bffile.params['bf']['range']):\n",
+    "    for krep, spiketimes in enumerate(bfspiketrains[bfs == stim]):\n",
+    "        bfspikecounts[kabi, krep] = within(spiketimes, bffile.params['delay'], bffile.params['delay'] + \\\n",
+    "                                        bffile.params['dur']*1)\n",
+    "#print(bfspikecounts.shape)\n",
+    "#fig, ax = plt.subplots(figsize = [5,4])\n",
+    "#\n",
+    "#ax.errorbar(bffile.params['bf']['range'], bfspikecounts.mean(axis=1), sps.stats.sem(bfspikecounts, axis=1))\n",
+    "#plt.xlabel('Frequency (Hz)')\n",
+    "#plt.ylabel('Spike Count')\n",
+    "##fig.savefig('ExampleFreq.png')\n",
+    "bfy = bfspikecounts.mean(axis=1)\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "peak_idx = np.argmax(bfy)\n",
+    "halfmax_y = (np.max(bfy) - np.min(bfy))/2 + np.min(bfy)\n",
+    "lup = np.where(bfy[:peak_idx] > halfmax_y)[0]\n",
+    "ldown = np.where(bfy[:peak_idx] < halfmax_y)[0]\n",
+    "rup = np.where(bfy[peak_idx+1:] > halfmax_y)[0]\n",
+    "rdown = np.where(bfy[peak_idx+1:] < halfmax_y)[0]\n",
+    "pm = np.array([-1,1])\n",
+    "\n",
+    "\n",
+    "if lup.size > 0 and ldown.size > 0:\n",
+    "    if between(halfmax_y, bfy[lup[0]] +  pm*bfse[lup[0]]):\n",
+    "        l = lup[1]\n",
+    "    elif between(halfmax_y, bfy[ldown[0]] +  pm*bfse[ldown[0]]):\n",
+    "        l = ldown[-1] \n",
+    "    elif lup[0] == ldown[-1] + 1:\n",
+    "        l = lup[0] \n",
+    "    else:\n",
+    "        l = lup[0]\n",
+    "elif lup.size > 0:\n",
+    "    l = lup[0]\n",
+    "elif ldown.size > 0: \n",
+    "    l = ldown[-1] \n",
+    "else:\n",
+    "    l = peak_idx\n",
+    "\n",
+    "if rup.size > 0 and rdown.size > 0:\n",
+    "    if between(halfmax_y, bfy[rup[-1] + peak_idx + 1] +  pm*bfse[rup[-1] + peak_idx + 1]):\n",
+    "        r = rup[-1] + peak_idx + 1\n",
+    "    elif between(halfmax_y, bfy[rdown[0] + peak_idx + 1] +  pm*bfse[rdown[0] + peak_idx + 1]):\n",
+    "        r = rdown[0] + peak_idx +1\n",
+    "    elif rup[-1] == rdown[0] - 1:\n",
+    "        r = rup[-1] + peak_idx + 1\n",
+    "    else:\n",
+    "        r = rup[-1] + peak_idx + 1\n",
+    "elif rup.size > 0:\n",
+    "    r = rup[-1] + peak_idx + 1\n",
+    "elif rdown.size > 0: \n",
+    "    r = rdown[0] + peak_idx +1\n",
+    "else:\n",
+    "    r = peak_idx   \n",
+    "\n",
+    "    \n",
+    "freq_lo, lofr = bfx[l], bfy[l]\n",
+    "freq_hi, hifr = bfx[r], bfy[r]\n",
+    "best_freq = np.mean([freq_lo, freq_hi])\n",
+    "display([best_freq, freq_lo, freq_hi])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#Best Freq\n",
+    "\n",
+    "fig, ax = plt.subplots(figsize = [5,4])\n",
+    "\n",
+    "ax.errorbar(bffile.params['bf']['range'], bfspikecounts.mean(axis=1), sps.stats.sem(bfspikecounts, axis=1))\n",
+    "plt.xlabel('Frequency (Hz)')\n",
+    "plt.ylabel('Spike Count')\n",
+    "\n",
+    "ax.scatter([freq_lo, freq_hi, best_freq], [lofr, hifr, np.max(bfy)], color = 'black')\n",
+    "plt.hlines((np.max(bfy) - np.min(bfy))/2 + np.min(bfy), 400, 10000)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#Frequency raster plot\n",
+    "fig, ax = plt.subplots()\n",
+    "ax.eventplot(bfspiketrains[bforder], color='k')\n",
+    "ax.set_yticks(np.arange(bffile.params['bf']['range'].size)[::4] \\\n",
+    "              * bffile.params['reps'] + bffile.params['reps']/2)\n",
+    "ax.set_yticklabels(bffile.params['bf']['range'][::4])\n",
+    "\n",
+    "ax.axvspan(bffile.params['delay'], bffile.params['delay']+bffile.params['dur']+10, color='yellow')\n",
+    "ax.axvspan(bffile.params['delay'], bffile.params['delay']+ 10, color='pink')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#ILD graph\n",
+    "fig, ax = plt.subplots(figsize= [5,4])\n",
+    "ax.errorbar(ildx, ildy, sps.stats.sem(ildspikecounts, axis=1), linewidth=2, color='black',marker = 'o')\n",
+    "ax.plot(ildx, Gaussian(ildx, *gopt), 'green')\n",
+    "ax.plot(ildx, Sigmoidal(ildx, *sopt), 'brown')\n",
+    "\n",
+    "plt.xlabel('ILD (dB)')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "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.9.12"
+  },
+  "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
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}

Shadron_Pena_Fig1.m

@@ -0,0 +1,503 @@
+%% Does the actual reliability calculations. Ruff removed HRTF's have less 
+% directions, so downsample normal files. Makes things go faster too.
+
+clear 
+%update these paths
+cd('F:\MATLAB\')
+addpath('F:\MATLAB\Brian_Reliability\Scripts')
+% % Load HRIRs
+
+owl_folder{1}='Merry';
+owl_folder{2}='Owl19';
+owl_folder{3}='Owl39';
+owl_folder{4}='Ugly';
+owl_folder{5}='Vespucci';
+
+haircut_id{1}='normal';
+haircut_id{2}='ohne'; % maximum cut
+date_folder = '06-06-22 Data'; %change to current date if rerunning next section
+%% Full analysis of reliability. Takes a while... 
+% Once this is done once, this and the next section can be skipped. Could
+% add another if loop to automatically skip this section.
+
+tic
+for o=2:5
+    for h=1:2
+        clearvars -except h haircut_id o owl_folder date_folder
+        load(['HRTF_Wagner_Lab\results\' owl_folder{o} '\' owl_folder{o} '__' haircut_id{h} '__hrir_result.mat'])
+        
+        newazimuths = -160:20:160;
+        newelevations = -60:20:80;
+        
+        dir=NaN(2,length(azimuths)*length(elevations));
+        TF1=NaN(length(azimuths)*length(elevations),size(hrir_l,3));
+        TF2=NaN(length(azimuths)*length(elevations),size(hrir_l,3));
+        count=0;
+        for a=1:length(azimuths)
+            for e=1:length(elevations)
+                count=count+1;
+                dir(:,count)=[elevations(e);azimuths(a)];
+                TF1(count,:)=hrir_l(a,e,:);
+                TF2(count,:)=hrir_r(a,e,:);
+            end
+        end        
+                
+        newdir=NaN(2,length(newazimuths)*length(newelevations));
+        count=0;
+        for a=1:length(newazimuths)
+            for e=1:length(newelevations)
+                count=count+1;
+                newdir(:,count)=[newelevations(e);newazimuths(a)];
+            end
+        end                
+        
+        dir = dir';
+        newdir = newdir';
+        downsample = ismember(dir, newdir, 'rows');
+        TF1 = TF1(downsample, :);
+        TF2 = TF2(downsample, :);
+        dir = newdir';
+        azimuth = newazimuths;
+        elevation = newelevations;
+        
+        clear newdir newazimuth newelevation
+        
+        % % Analysis parameters
+        
+        %center frequencies of cochlear filters (Hz)
+        cF = 400:200:10000;
+        %cF = linspace(1000,9000, 9);
+        n_cF=length(cF);
+        
+        %Sampling frequency for HRTFs (Hz)
+        Fs = samplingrate;
+        
+        %Factor by which you will upsample
+        Factor = 5;
+        %Number of trials for each condition
+        Nt=5;
+        % % Get cochlear filters
+        %Upsampled frequency
+        FS=Fs*Factor;
+        %Get filters
+        fcoefs=getFilterCoefs(cF,FS);
+        
+        
+        % % Number of directions and initialize        
+        nd_tot=size(TF1,1);
+        n_dir = size(TF1,1);
+        
+        IPDm_all = zeros(n_dir, n_cF, Nt);
+        IPDs_all = zeros(n_dir, n_cF, Nt);
+%         ITDm=zeros(n_dir,n_cF);
+%         ITDs=zeros(n_dir,n_cF);
+        ILDm_all = zeros(n_dir, n_cF, Nt);
+        ILDs_all = zeros(n_dir, n_cF, Nt);
+        
+        gainr_all = zeros(n_dir, n_cF, Nt);
+        gainl_all = zeros(n_dir, n_cF, Nt);
+        
+        
+        % %
+        tic
+        
+        TF1_res = resample(TF1',Factor,1)';
+        TF2_res = resample(TF2',Factor,1)';
+        
+        parfor idir = 1:n_dir % TODO parfor
+            
+            %Resample the HRIRs
+            hL=TF1_res(idir,:);
+            hR=TF2_res(idir,:);
+            %Run over trials
+            
+            for tt=1:Nt
+                
+                %Generate a target sound
+                s=genSignal(100,1/30,1,2,2*pi*12);
+                s1=resample(s,Factor,1);
+                
+                %Generate a masker
+                s=genSignal(100,1/30,1,2,2*pi*12);
+                s2=resample(s,Factor,1);
+                
+                %Convolve target with HRIR
+                sL1=conv(hL,s1);
+                sR1=conv(hR,s1);
+                
+                Z=zeros(n_dir,n_cF);
+                GL=zeros(n_dir,n_cF);
+                GR=zeros(n_dir,n_cF);
+                ITDp=zeros(n_cF,n_dir);
+                
+                %Run over masking directions
+                %disp(num2str([toc o h idir tt]))
+                
+                for k=1:nd_tot
+                    %disp(num2str([o h idir tt k]))
+                    
+                    %Resample the HRIRs
+                    hL2=TF1_res(k,:); %#ok<PFBNS>
+                    hR2=TF2_res(k,:); %#ok<PFBNS>
+                    
+                    %Convolve masker with HRIR and add to target
+                    sL=sL1+conv(hL2,s2);
+                    sR=sR1+conv(hR2,s2);
+                    
+                    %cochlear (gammmatone) filterbank
+                    vL=ERBFilterBank(sL,fcoefs);
+                    vR=ERBFilterBank(sR,fcoefs);
+                    
+                    %Run over frequency and compute cross correlation
+                    for icF = 1:n_cF
+                        [x,l]=xcorr(vL(icF,:),vR(icF,:));
+                        [~,j]=max(x);
+                        ITDp(icF,k)=l(j)*1/Fs*1e6/Factor;
+                    end
+                    
+                    %Compute ILD
+                    R=rms(vR,2);
+                    L=rms(vL,2);
+                    Z(k,:)=20*log10(R./L);
+                    GR(k,:)=20*log10(R);
+                    GL(k,:)=20*log10(L);
+                    
+                end
+                
+                
+                %Run over frequency and compute mean and SD
+                for n=1:n_cF
+                 
+                    IPDp=ITDp(n,:)*cF(n)/1e6*2*pi; %#ok<PFBNS>
+                    IPDs_all(idir,n,tt) = circstd(IPDp/Nt);
+                    IPDm_all(idir,n,tt) = circmean(IPDp/Nt);
+                    
+                end
+                
+                %compute mean and SD for ILD. 
+                ILDm_all(idir,:,tt) = mean(Z,1);
+                ILDs_all(idir,:,tt) = std(Z,[],1);
+                
+                %compute gain
+                gainr_all(idir,:,tt) = mean(GR,1);
+                gainl_all(idir,:,tt) = mean(GL,1);
+
+                
+            end     
+        end
+        IPDs = sum(IPDs_all, 3);
+        IPDm = sum(IPDm_all, 3);
+        ITDs = (1e6 * IPDs) ./ (2 * pi * cF);
+        ITDm = (1e6 * IPDm) ./ (2 * pi * cF);
+        ILDs = sum(ILDs_all, 3)/Nt;
+        ILDm = sum(ILDm_all, 3)/Nt;
+        gainl = sum(gainl_all, 3)/Nt;
+        gainr = sum(gainr_all, 3)/Nt;
+
+        toc
+        save(['Shadron_Pena_2022\', date_folder, '\', owl_folder{o} '__' haircut_id{h} '__hrtf_fig1.mat'])
+    end
+end
+%% Some minor editing
+
+for o = 1:5
+    for h = 1:2
+        clearvars -except h haircut_id o owl_folder date_folder
+        load(['Shadron_Pena_2022\', date_folder, '\', owl_folder{o} '__' haircut_id{h} '__hrtf_fig1.mat'])        
+        
+        %have to wrap ITD's to make things make sense
+        %Earlier had to convert 17x8x563 to a 136x563, so now have to revert
+        %back to 17x8x20.
+        
+        x = size(hrir_l, 3);
+        y = size(cF, 2);
+            
+            IxD.IPDm = zeros(17,8,y);
+            IxD.IPDs = zeros(17,8,y);
+            IxD.ITDm = zeros(17,8,y);
+            IxD.ITDs = zeros(17,8,y);
+            IxD.ILDm = zeros(17,8,y);
+            IxD.ILDs = zeros(17,8,y);
+            IxD.gainl = zeros(17,8,y);
+            IxD.gainr = zeros(17,8,y);
+            
+            for n = 1:y
+                for a = 1:17
+                    for e = 1:8
+                        counter = (a-1)*8 + e;
+                        IxD.IPDm(a,e,n) = IPDm(counter,n);
+                        IxD.IPDs(a,e,n) = IPDs(counter,n);
+                        IxD.ITDm(a,e,n) = ITDm(counter,n);
+                        IxD.ITDs(a,e,n) = ITDs(counter,n);
+                        IxD.ILDm(a,e,n) = ILDm(counter,n);
+                        IxD.ILDs(a,e,n) = ILDs(counter,n);
+                        IxD.gainl(a,e,n) = gainl(counter,n);
+                        IxD.gainr(a,e,n) = gainr(counter,n);
+                    end
+                end
+            end
+
+        clear n a e count Factor Fs FS idir n n_cF n_dir nd_tot Nt tt x ILDm ILDs IPDm IPDp IPDs ITDp ITDs ITDm
+        clear samplingrate s1 s2 sL1 sR1 TF1 TF2 Z ax counter fcoef dir s fcoefs hL hR hrir_l hrir_r 
+        
+        if size(IxD.IPDm,2) > 8 %Loop occurs if spatial locations are > 17 azimuth and > 9 elevation positions
+            %Resamples to have locations that are
+            %-160:20:160 in azimuth and -60:20:80 in elevation
+            for n = 1:y
+                for a = 1:2:34
+                    for e = 2:2:16
+                        IxD2.IPDm((a+1)/2, e/2, n) = IxD.IPDm(a,e,n);
+                        IxD2.IPDs((a+1)/2, e/2, n) = IxD.IPDs(a,e,n);
+                        IxD2.ITDm((a+1)/2, e/2, n) = IxD.ITDm(a,e,n);
+                        IxD2.ITDs((a+1)/2, e/2, n) = IxD.ITDs(a,e,n);
+                        IxD2.ILDm((a+1)/2, e/2, n) = IxD.ILDm(a,e,n);
+                        IxD2.ILDs((a+1)/2, e/2, n) = IxD.ILDs(a,e,n);
+                    end
+                end
+            end
+            
+            
+            clear IxD
+            IxD = IxD2;
+            clear IxD2
+            
+            azimuths = -160:20:160;
+            elevations = -60:20:80;
+        end
+        
+        
+        save(['Shadron_Pena_2022\', date_folder, '\', owl_folder{o} '__' haircut_id{h} '__stats_fig1.mat'])
+        
+    end
+end
+%% Compile all data into one struct
+clearvars -except h haircut_id o owl_folder date_folder
+
+condition = [];
+for o = [1 2 3 4 5]
+    for h = 1:2
+        clearvars -except h haircut_id o owl_folder condition date_folder
+        load(['Shadron_Pena_2022\', date_folder, '\', owl_folder{o} '__' haircut_id{h} '__stats_fig1.mat'])
+        condition.(owl_folder{o}).(haircut_id{h}) = IxD;
+         
+    end
+end
+            
+clear IxD name
+
+%average the owls
+for h=1:2
+for o = [1 2 3 4 5]
+        if o == 1
+            condition.(haircut_id{h}).avg.IPDm = condition.(owl_folder{o}).(haircut_id{h}).IPDm;
+            condition.(haircut_id{h}).avg.IPDs = condition.(owl_folder{o}).(haircut_id{h}).IPDs;
+            condition.(haircut_id{h}).avg.ITDm = condition.(owl_folder{o}).(haircut_id{h}).ITDm;
+            condition.(haircut_id{h}).avg.ITDs = condition.(owl_folder{o}).(haircut_id{h}).ITDs;
+            condition.(haircut_id{h}).avg.ILDm = condition.(owl_folder{o}).(haircut_id{h}).ILDm;
+            condition.(haircut_id{h}).avg.ILDs = condition.(owl_folder{o}).(haircut_id{h}).ILDs;
+            condition.(haircut_id{h}).avg.gainl = condition.(owl_folder{o}).(haircut_id{h}).gainl;
+            condition.(haircut_id{h}).avg.gainr = condition.(owl_folder{o}).(haircut_id{h}).gainr;
+
+        else
+            condition.(haircut_id{h}).avg.IPDs = condition.(haircut_id{h}).avg.IPDs + condition.(owl_folder{o}).(haircut_id{h}).IPDs;
+            condition.(haircut_id{h}).avg.ITDm = condition.(haircut_id{h}).avg.ITDm + condition.(owl_folder{o}).(haircut_id{h}).ITDm;
+            condition.(haircut_id{h}).avg.ITDs = condition.(haircut_id{h}).avg.ITDs + condition.(owl_folder{o}).(haircut_id{h}).ITDs;
+            condition.(haircut_id{h}).avg.ILDm = condition.(haircut_id{h}).avg.ILDm + condition.(owl_folder{o}).(haircut_id{h}).ILDm;
+            condition.(haircut_id{h}).avg.ILDs = condition.(haircut_id{h}).avg.ILDs + condition.(owl_folder{o}).(haircut_id{h}).ILDs;
+            condition.(haircut_id{h}).avg.gainl = condition.(haircut_id{h}).avg.gainl + condition.(owl_folder{o}).(haircut_id{h}).gainl;
+            condition.(haircut_id{h}).avg.gainr = condition.(haircut_id{h}).avg.gainr + condition.(owl_folder{o}).(haircut_id{h}).gainr;
+        end
+end
+condition.(haircut_id{h}).avg.IPDm = condition.(haircut_id{h}).avg.IPDm / 5;
+condition.(haircut_id{h}).avg.IPDs = condition.(haircut_id{h}).avg.IPDs / 5;
+condition.(haircut_id{h}).avg.ITDm = condition.(haircut_id{h}).avg.ITDm / 5;
+condition.(haircut_id{h}).avg.ITDs = condition.(haircut_id{h}).avg.ITDs / 5;
+condition.(haircut_id{h}).avg.ILDm = condition.(haircut_id{h}).avg.ILDm / 5;
+condition.(haircut_id{h}).avg.ILDs = condition.(haircut_id{h}).avg.ILDs / 5;
+condition.(haircut_id{h}).avg.gainl = condition.(haircut_id{h}).avg.gainl/5;
+condition.(haircut_id{h}).avg.gainr = condition.(haircut_id{h}).avg.gainr/5;
+end
+
+%% Makes the reliability maps across azimuth as in Cazettes et al (2014)
+
+lowf = 1200; % The lowest frequency you want to include in the normalization
+highf = 8000;
+%low = ((lowf-1000)/200); 
+low = find(cF == lowf) - 1;
+high = find(cF == highf);
+
+for h = 1:2
+    condition.(haircut_id{h}).avg.reliability.IPD = (condition.(haircut_id{h}).avg.IPDs.^-1);
+    condition.(haircut_id{h}).avg.reliability.ILD = (condition.(haircut_id{h}).avg.ILDs.^-1);
+    condition.(haircut_id{h}).avg.reliability.normalizedIPD = zeros(size(azimuths,1),high-low);
+    condition.(haircut_id{h}).avg.reliability.normalizedILD = zeros(size(azimuths,1),high-low);
+    condition.(haircut_id{h}).avg.reliability.normIPDaz = zeros(size(azimuths,1),2);
+
+    
+    for a = 1:17
+        [big, pt] = max(condition.(haircut_id{h}).avg.reliability.IPD(a,4,1+low:high));
+        for n = 1:high-low
+            condition.(haircut_id{h}).avg.reliability.normalizedIPD(a,n) = condition.(haircut_id{h}).avg.reliability.IPD(a,4,n+low)/big;
+        end
+        condition.(haircut_id{h}).avg.reliability.normIPDaz(a,1) = mean(condition.(haircut_id{h}).avg.reliability.IPD(a,4,1+low:high).^-1);
+        condition.(haircut_id{h}).avg.reliability.normIPDaz(a,2) = std(condition.(haircut_id{h}).avg.reliability.IPD(a,4,1+low:high).^-1);
+    end
+    for a = 1:17
+        [big, pt] = max(condition.(haircut_id{h}).avg.reliability.ILD(a,4,1+low:high));
+        for n = 1:high-low
+            condition.(haircut_id{h}).avg.reliability.normalizedILD(a,n) = condition.(haircut_id{h}).avg.reliability.ILD(a,4,n+low)/big;
+        end
+    end
+
+end
+
+short = cF(1+low:high);
+
+% Can run to see overall pattern before normalization
+% figure
+% normalmap = imagesc(azimuths, short, condition.(haircut_id{1}).avg.reliability.normalizedIPD'); axis xy
+% title(['IPD Reliability Along Azimuth: ', haircut_id{1}, ' (n = 5)']); colorbar; colormap jet
+% 
+% figure
+% ohnemap = imagesc(azimuths, short, condition.(haircut_id{2}).avg.reliability.normalizedIPD'); axis xy
+% title(['IPD Reliability Along Azimuth: ', haircut_id{2}, ' (n = 5)']); colorbar; colormap jet
+
+%% Figure 1a and 1b
+
+posazi = 0:20:160;
+shortKHz = short/1000;
+
+colors = f_gsafecmap(256);
+colormap(flipud(colors));
+
+d = figure;
+set(gcf,'Position',[100 100 600 450]);
+posnormalmap = imagesc(posazi, shortKHz, condition.(haircut_id{1}).avg.reliability.normalizedIPD(9:17,:)'); axis xy
+title('Normal')
+%title(['IPD Reliability Along Azimuth: ', haircut_id{1}, ' (n = 5)']); 
+posnormalmapc = colorbar; colormap(flipud(colors)); clim([0.2 1]);
+%posnormalmapc.Label.String = 'Reliability'; posnormalmapc.Label.FontSize = 18;
+ylabel('Frequency (KHz)'); xlabel('Azimuth (°)');
+ylabel('Frequency (KHz)', 'FontSize', 22); xlabel('Azimuth (°)', 'FontSize', 22);
+set(gca,'XTick',[0 45 90 135 180], 'XTickLabel', [0 45 90 135 180], 'fontsize', 18);
+set(d,'Units','Inches');
+set(d,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[2.8, 2.24])
+text(-35,8.6, 'a', 'FontSize', 28);
+%saveas(posnormalmap, 'F:\Various Files\My Papers\Shadron and Pena 2022\Matlab\normreliab.png')
+
+figure
+set(gcf,'Position',[100 100 600 450]);
+posohnemap = imagesc(posazi, shortKHz, condition.(haircut_id{2}).avg.reliability.normalizedIPD(9:17,:)'); axis xy
+title('Ruffcut')
+%title(['IPD Reliability Along Azimuth: ', haircut_id{2}, ' (n = 5)']); 
+posohnemapc = colorbar; colormap(flipud(colors)); clim([0.2 1]);
+posohnemapc.Label.String = 'Reliability'; posohnemapc.Label.FontSize = 22;
+%ylabel('Frequency (KHz)'); xlabel('Azimuth (deg)');
+xlabel('Azimuth (°)', 'FontSize', 22);
+set(gca,'XTick',[0 45 90 135 180], 'XTickLabel', [0 45 90 135 180], 'fontsize', 18);
+text(-35,8.6, 'b', 'FontSize', 28);
+%saveas(posohnemap, 'F:\Various Files\My Papers\Shadron and Pena 2022\Matlab\ruffcutreliab.png')
+
+
+%For non normalized:
+%imagesc(posazi, shortKHz, permute(squeeze(condition.(haircut_id{1}).avg.reliability.IPD(9:17,4,:)),[2 1]));colorbar
+%% Figure 1c
+
+diffIPD = condition.(haircut_id{2}).avg.reliability.IPD - condition.(haircut_id{1}).avg.reliability.IPD;
+diffILD = condition.(haircut_id{2}).avg.reliability.ILD - condition.(haircut_id{1}).avg.reliability.ILD;
+dIPD = squeeze(diffIPD(:,4,:));
+dILD = squeeze(diffILD(:,4,:));
+
+figure
+set(gcf,'Position',[100 100 600 450]);
+IPDd = imagesc(azimuths(9:end), cF(low:high), dIPD(9:end,low:high)'); axis xy; 
+IPDdc = colorbar; clim([-.3 .3]); colors = f_gsafecmap(256); colormap(flipud(colors));
+ylabel('Frequency (KHz)', 'FontSize', 22); 
+xlabel('Azimuth (°)', 'FontSize', 22);
+IPDdc.Label.String = 'Ruffcut IPD - Normal IPD';
+IPDdc.Label.FontSize = 22;
+IPDdc.FontSize = 16;
+set(gca,'XTick',[0 45 90 135 180], 'XTickLabel', [0 45 90 135 180], 'fontsize', 18);
+set(gca,'YTick',2000:2000:8000, 'YTickLabel', [2 4 6 8], 'fontsize', 18);
+text(-35,8600, 'c', 'FontSize', 24);
+%saveas(IPDd, 'F:\Various Files\My Papers\Shadron and Pena 2022\Matlab\reliabdiff.png')
+
+%set(IPDdc,'Units','Inches');
+
+% colors = f_gsafecmap(256);
+% colormap(flipud(colors));
+% ILDd = imagesc(azimuths(9:end), cF(4:end), dILD(9:end,4:end)'); axis xy; colorbar; caxis([-.5 .5]);
+% ylabel('Frequency (KHz)','FontSize', 20); 
+% xlabel('Azimuth (Deg)','FontSize', 20);
+% title('Difference in ILD reliability', 'FontSize', 20)
+
+%% ILD, not shown in manuscript
+
+posazi = 160:-20:0;
+shortKHz = short/1000;
+
+colors = f_gsafecmap(256);
+colormap(flipud(colors));
+
+figure
+negnormalmap = imagesc(posazi, shortKHz, condition.(haircut_id{1}).avg.reliability.normalizedILD(1:9,:)'); axis xy
+title('Normal', 'FontSize', 20); 
+negnormalmapc = colorbar;  clim([0.2 1]);
+negnormalmapc.Label.String = 'Reliability'; negnormalmapc.Label.FontSize = 18;
+colormap(flipud(colors));
+ylabel('Frequency (KHz)', 'FontSize', 20); xlabel('Azimuth (deg)', 'FontSize', 20);
+set(gca,'XTick',[0 45 90 135 180], 'XTickLabel', [0 45 90 135 180], 'fontsize', 16);
+
+figure
+negohnemap = imagesc(posazi, shortKHz, condition.(haircut_id{2}).avg.reliability.normalizedILD(1:9,:)'); axis xy
+title('Ruff Cut', 'FontSize', 20); 
+negohnemapc = colorbar; colormap(flipud(colors)); clim([0.2 1]);
+negohnemapc.Label.String = 'Reliability'; negohnemapc.Label.FontSize = 18;
+%ylabel('Frequency (KHz)','FontSize', 20); 
+xlabel('Azimuth (deg)', 'FontSize', 20);
+set(gca,'XTick',[0 45 90 135 180], 'XTickLabel', [0 45 90 135 180], 'fontsize', 16);
+
+%For non normalized:
+%imagesc(posazi, shortKHz, permute(squeeze(condition.(haircut_id{1}).avg.reliability.IPD(9:17,4,:)),[2 1]));colorbar
+
+%% ILD reliability, not used in manuscript
+colors = f_gsafecmap(256);
+colormap(flipud(colors));
+
+figure
+posnormalmap = imagesc(posazi, shortKHz, condition.(haircut_id{1}).avg.reliability.normalizedILD(9:17,:)'); axis xy
+title('Normal: ILD reliability')
+%title(['ILD Reliability Along Azimuth: ', haircut_id{1}, ' (n = 5)']); 
+posnormalmapc = colorbar; colormap(flipud(colors));
+posnormalmapc.Label.String = 'Reliability'; posnormalmapc.Label.FontSize = 12;
+ylabel('Frequency (KHz)'); xlabel('ITD (us)');
+
+figure
+posohnemap = imagesc(posazi, shortKHz, condition.(haircut_id{2}).avg.reliability.normalizedILD(9:17,:)'); axis xy
+title('Ruffcut: ILD reliability')
+%title(['IPD Reliability Along Azimuth: ', haircut_id{2}, ' (n = 5)']); 
+posohnemapc = colorbar; colormap(flipud(colors));
+posohnemapc.Label.String = 'Reliability'; posohnemapc.Label.FontSize = 12;
+ylabel('Frequency (KHz)'); xlabel('ITD (us)');
+
+%For non normalized:
+%imagesc(posazi, shortKHz, permute(squeeze(condition.(haircut_id{1}).avg.reliability.IPD(9:17,4,:)),[2 1]));colorbar
+
+%% s.d. IPD averaged across freq for each azimuth, not used in manuscript
+
+IPDstd = figure;
+hold on
+scatter(azimuths, condition.(haircut_id{1}).avg.reliability.normIPDaz(:,1)/180, 40, 'blue', 'filled', 'o')
+scatter(azimuths, condition.(haircut_id{2}).avg.reliability.normIPDaz(:,1)/180, 40, 'red', 'filled', 'o')
+
+
+errorbar(azimuths, condition.(haircut_id{1}).avg.reliability.normIPDaz(:,1)/180, condition.(haircut_id{1}).avg.reliability.normIPDaz(:,2)/180, ...
+     'blue', 'LineStyle', 'none')
+errorbar(azimuths, condition.(haircut_id{2}).avg.reliability.normIPDaz(:,1)/180, condition.(haircut_id{2}).avg.reliability.normIPDaz(:,2)/180, ...
+     'red', 'LineStyle', 'none')
+xlabel('Azimuth (deg)')
+ylabel('IPD standard deviation')
+xlim([-10,90])
+
+legend('Normal', 'Ruff-Removed', 'Location', 'southeast')
+fontsize(gca, 16,"points")

Shadron_Pena_Fig5.m

@@ -0,0 +1,460 @@
+%% Does the actual reliability calculations. Ruff removed HRTF's have less 
+% directions, so downsample normal files. Makes things go faster too.
+% Currently at 10 iterations
+clear 
+cd('F:\Shadron_Pena_2022\')
+%cd('C:\Users\penalab2\Documents\MATLAB\Brian_Reliability')
+addpath('F:\Shadron_Pena_2022')
+% % Load HRIRs
+
+owl_folder{1}='Merry';
+owl_folder{2}='Owl19';
+owl_folder{3}='Owl39';
+owl_folder{4}='Ugly';
+owl_folder{5}='Vespucci';
+
+haircut_id{1}='normal';
+haircut_id{2}='ohne'; % maximum cut
+date_folder = '10-24-22 Data'; %change to current date
+%% This anlaysis is similar to Figure 1 analysis, but no maskers are used.
+
+for o=1:5
+    for h=1:2
+        clearvars -except h haircut_id o owl_folder date_folder
+        load(['HRTF_Wagner_Lab\results\' owl_folder{o} '\' owl_folder{o} '__' haircut_id{h} '__hrir_result.mat'])
+        
+        newazimuths = -160:20:160;
+        newelevations = -60:20:80;
+        
+        dir=NaN(2,length(azimuths)*length(elevations));
+        TF1=NaN(length(azimuths)*length(elevations),size(hrir_l,3));
+        TF2=NaN(length(azimuths)*length(elevations),size(hrir_l,3));
+        count=0;
+        for a=1:length(azimuths)
+            for e=1:length(elevations)
+                count=count+1;
+                dir(:,count)=[elevations(e);azimuths(a)];
+                TF1(count,:)=hrir_l(a,e,:);
+                TF2(count,:)=hrir_r(a,e,:);
+            end
+        end        
+                
+        newdir=NaN(2,length(newazimuths)*length(newelevations));
+        count=0;
+        for a=1:length(newazimuths)
+            for e=1:length(newelevations)
+                count=count+1;
+                newdir(:,count)=[newelevations(e);newazimuths(a)];
+            end
+        end                
+        
+        dir = dir';
+        newdir = newdir';
+        downsample = ismember(dir, newdir, 'rows');
+        TF1 = TF1(downsample, :);
+        TF2 = TF2(downsample, :);
+        dir = newdir';
+        azimuth = newazimuths;
+        elevation = newelevations;
+        
+        clear newdir newazimuth newelevation
+        
+        % % Analysis parameters
+        
+        %center frequencies of cochlear filters (Hz)
+        cF = 400:200:10000;
+        %cF = linspace(1000,9000, 9);
+        n_cF=length(cF);
+        
+        %Sampling frequency for HRTFs (Hz)
+        Fs = samplingrate;
+        
+        %Factor by which you will upsample
+        Factor = 5;
+        %Number of trials for each condition
+        Nt=10;
+        % % Get cochlear filters
+        %Upsampled frequency
+        FS=Fs*Factor;
+        %Get filters
+        fcoefs=getFilterCoefs(cF,FS);
+
+
+        % % Number of directions and initialize
+        n_dir = size(TF1,1);
+
+        IPDm_all = zeros(n_dir, n_cF, Nt);
+        IPDs_all = zeros(n_dir, n_cF, Nt);
+        %         ITDm=zeros(n_dir,n_cF);
+        %         ITDs=zeros(n_dir,n_cF);
+        ILDm_all = zeros(n_dir, n_cF, Nt);
+        ILDs_all = zeros(n_dir, n_cF, Nt);
+
+        gainr_all = zeros(n_dir, n_cF, Nt);
+        gainl_all = zeros(n_dir, n_cF, Nt);
+
+
+        % %
+        tic
+
+        TF1_res = resample(TF1',Factor,1)';
+        TF2_res = resample(TF2',Factor,1)';
+
+        parfor idir = 1:n_dir % TODO parfor
+
+            %Resample the HRIRs
+            hL=TF1_res(idir,:);
+            hR=TF2_res(idir,:);
+            %Run over trials
+
+            for tt=1:Nt
+
+                %Generate a target sound
+                s=genSignal(100,1/30,1,2,2*pi*12);
+                s1=resample(s,Factor,1);
+
+                %Convolve target with HRIR
+                sL=conv(hL,s1);
+                sR=conv(hR,s1);
+
+                ITDp=zeros(n_cF);
+
+                %cochlear (gammmatone) filterbank
+                vL=ERBFilterBank(sL,fcoefs);
+                vR=ERBFilterBank(sR,fcoefs);
+
+                %Run over frequency and compute cross correlation
+                for icF = 1:n_cF
+                    [itd,l]=xcorr(vL(icF,:),vR(icF,:));
+                    [~,j]=max(itd);
+                    ITDp(icF)=l(j)*1/Fs*1e6/Factor;
+                end
+
+                %Compute ILD
+                R=rms(vR,2);
+                L=rms(vL,2);
+                Z=20*log10(R./L);
+                GR=20*log10(R);
+                GL=20*log10(L);
+                
+                %Run over frequency and compute mean and SD
+                for n=1:n_cF
+                    IPD_all(idir,n,tt) =ITDp(n)*cF(n)/1e6*2*pi; %#ok<PFBNS>
+                end
+
+                %compute mean and SD for ILD.
+                ILD_all(idir,:,tt) = Z;
+
+                %compute gain
+                gainr_all(idir,:,tt) =GR;
+                gainl_all(idir,:,tt) =GL;
+
+
+            end
+        end
+        IPDs = std(IPD_all,0,3);
+        IPDm = mean(IPD_all, 3);
+        ITDs = (1e6 * IPDs) ./ (2 * pi * cF);
+        ITDm = (1e6 * IPDm) ./ (2 * pi * cF);
+        ILDs = std(ILD_all,0,3);
+        ILDm = sum(ILD_all, 3)/Nt;
+        gainl = sum(gainl_all, 3)/Nt;
+        gainr = sum(gainr_all, 3)/Nt;
+
+        toc
+        save(['Shadron_Pena_2022\', date_folder, '\', owl_folder{o} '__' haircut_id{h} '__hrtf_fig5.mat'])
+    end
+end
+%% Some minor editing
+
+for o = 1:5
+    for h = 1:2
+        clearvars -except h haircut_id o owl_folder date_folder
+        load(['Shadron_Pena_2022\', date_folder, '\', owl_folder{o} '__' haircut_id{h} '__hrtf_fig5.mat'])        
+        
+        %have to wrap ITD's to make things make sense
+        %Earlier had to convert 17x8x563 to a 136x563, so now have to revert
+        %back to 17x8x20.
+        
+        itd = size(hrir_l, 3);
+        y = size(cF, 2);
+            
+            IxD.IPDm = zeros(17,8,y);
+            IxD.IPDs = zeros(17,8,y);
+            IxD.ITDm = zeros(17,8,y);
+            IxD.ITDs = zeros(17,8,y);
+            IxD.ILDm = zeros(17,8,y);
+            IxD.ILDs = zeros(17,8,y);
+            IxD.gainl = zeros(17,8,y);
+            IxD.gainr = zeros(17,8,y);
+            
+            for n = 1:y
+                for a = 1:17
+                    for e = 1:8
+                        counter = (a-1)*8 + e;
+                        IxD.IPDm(a,e,n) = IPDm(counter,n);
+                        IxD.IPDs(a,e,n) = IPDs(counter,n);
+                        IxD.ITDm(a,e,n) = ITDm(counter,n);
+                        IxD.ITDs(a,e,n) = ITDs(counter,n);
+                        IxD.ILDm(a,e,n) = ILDm(counter,n);
+                        IxD.ILDs(a,e,n) = ILDs(counter,n);
+                        IxD.gainl(a,e,n) = gainl(counter,n);
+                        IxD.gainr(a,e,n) = gainr(counter,n);
+                    end
+                end
+            end
+
+        clear n a e count Factor Fs FS idir n n_cF n_dir Nt tt itd ILDm ILDs IPDm IPDp IPDs ITDp ITDs ITDm
+        clear samplingrate s1 s2 sL1 sR1 TF1 TF2 Z ax counter fcoef dir s fcoefs hL hR hrir_l hrir_r 
+        
+        if size(IxD.IPDm,2) > 8 %Loop occurs if spatial locations are > 17 azimuth and > 9 elevation positions
+            %Resamples to have locations that are
+            %-160:20:160 in azimuth and -60:20:80 in elevation
+            for n = 1:y
+                for a = 1:2:34
+                    for e = 2:2:16
+                        IxD2.IPDm((a+1)/2, e/2, n) = IxD.IPDm(a,e,n);
+                        IxD2.IPDs((a+1)/2, e/2, n) = IxD.IPDs(a,e,n);
+                        IxD2.ITDm((a+1)/2, e/2, n) = IxD.ITDm(a,e,n);
+                        IxD2.ITDs((a+1)/2, e/2, n) = IxD.ITDs(a,e,n);
+                        IxD2.ILDm((a+1)/2, e/2, n) = IxD.ILDm(a,e,n);
+                        IxD2.ILDs((a+1)/2, e/2, n) = IxD.ILDs(a,e,n);
+                    end
+                end
+            end
+            
+            
+            clear IxD
+            IxD = IxD2;
+            clear IxD2
+            
+            azimuths = -160:20:160;
+            elevations = -60:20:80;
+        end
+        
+        
+        save(['Shadron_Pena_2022\', date_folder, '\', owl_folder{o} '__' haircut_id{h} '__stats_fig5.mat'])
+        
+    end
+end
+%% Compile all data into one struct
+clearvars -except h haircut_id o owl_folder date_folder
+
+condition = [];
+for o = [1 2 3 4 5]
+    for h = 1:2
+        clearvars -except h haircut_id o owl_folder condition date_folder
+        load(['Shadron_Pena_2022\', date_folder, '\', owl_folder{o} '__' haircut_id{h} '__stats_fig5.mat'])
+        condition.(owl_folder{o}).(haircut_id{h}) = IxD;
+         
+    end
+end
+            
+clear IxD name
+
+%average the owls
+for h=1:2
+for o = [1 2 3 4 5]
+        if o == 1
+            condition.(haircut_id{h}).avg.IPDm = condition.(owl_folder{o}).(haircut_id{h}).IPDm;
+            condition.(haircut_id{h}).avg.IPDs = condition.(owl_folder{o}).(haircut_id{h}).IPDs;
+            condition.(haircut_id{h}).avg.ITDm = condition.(owl_folder{o}).(haircut_id{h}).ITDm;
+            condition.(haircut_id{h}).avg.ITDs = condition.(owl_folder{o}).(haircut_id{h}).ITDs;
+            condition.(haircut_id{h}).avg.ILDm = condition.(owl_folder{o}).(haircut_id{h}).ILDm;
+            condition.(haircut_id{h}).avg.ILDs = condition.(owl_folder{o}).(haircut_id{h}).ILDs;
+            condition.(haircut_id{h}).avg.gainl = condition.(owl_folder{o}).(haircut_id{h}).gainl;
+            condition.(haircut_id{h}).avg.gainr = condition.(owl_folder{o}).(haircut_id{h}).gainr;
+
+        else
+            condition.(haircut_id{h}).avg.IPDs = condition.(haircut_id{h}).avg.IPDs + condition.(owl_folder{o}).(haircut_id{h}).IPDs;
+            condition.(haircut_id{h}).avg.ITDm = condition.(haircut_id{h}).avg.ITDm + condition.(owl_folder{o}).(haircut_id{h}).ITDm;
+            condition.(haircut_id{h}).avg.ITDs = condition.(haircut_id{h}).avg.ITDs + condition.(owl_folder{o}).(haircut_id{h}).ITDs;
+            condition.(haircut_id{h}).avg.ILDm = condition.(haircut_id{h}).avg.ILDm + condition.(owl_folder{o}).(haircut_id{h}).ILDm;
+            condition.(haircut_id{h}).avg.ILDs = condition.(haircut_id{h}).avg.ILDs + condition.(owl_folder{o}).(haircut_id{h}).ILDs;
+            condition.(haircut_id{h}).avg.gainl = condition.(haircut_id{h}).avg.gainl + condition.(owl_folder{o}).(haircut_id{h}).gainl;
+            condition.(haircut_id{h}).avg.gainr = condition.(haircut_id{h}).avg.gainr + condition.(owl_folder{o}).(haircut_id{h}).gainr;
+        end
+end
+condition.(haircut_id{h}).avg.IPDm = condition.(haircut_id{h}).avg.IPDm / 5;
+condition.(haircut_id{h}).avg.IPDs = condition.(haircut_id{h}).avg.IPDs / 5;
+condition.(haircut_id{h}).avg.ITDm = condition.(haircut_id{h}).avg.ITDm / 5;
+condition.(haircut_id{h}).avg.ITDs = condition.(haircut_id{h}).avg.ITDs / 5;
+condition.(haircut_id{h}).avg.ILDm = condition.(haircut_id{h}).avg.ILDm / 5;
+condition.(haircut_id{h}).avg.ILDs = condition.(haircut_id{h}).avg.ILDs / 5;
+condition.(haircut_id{h}).avg.gainl = condition.(haircut_id{h}).avg.gainl/5;
+condition.(haircut_id{h}).avg.gainr = condition.(haircut_id{h}).avg.gainr/5;
+end
+
+
+
+%% Code for constructing Figure 5, model neurons
+
+%model neuron input parameters 
+freq = 2000:1000:7000;
+width = length(freq);
+sigma = length(freq);
+itd = -800:800;
+ild = -20:0.01:20;
+
+%tuning for model neurons
+itd_peak = 0;
+ild_peak = 0;
+
+%max spike count 
+A = 10;
+
+cF_ind = 9:5:34;
+
+normalITD = zeros(length(azimuths), length(elevations), length(freq));
+normalILD = zeros(length(azimuths), length(elevations), length(freq));
+ruffcutITD = zeros(length(azimuths), length(elevations), length(freq));
+ruffcutILD = zeros(length(azimuths), length(elevations), length(freq));
+
+for p = 1:length(freq)
+    
+    %model ITD curve
+    itd_y = A* (exp(cos( (2*pi*freq(p)/1000000) *(itd - itd_peak))) - exp(-1)      ) / (exp(1) - exp(-1));
+    
+    %model ILD curve
+    ild_y = 10* exp( -(8)^-1* (ild - ild_peak).^2);
+
+    %determine spike count given the ITD from HRTFs
+    for q = 1:numel(condition.normal.avg.ITDm(:,:,p))
+        [row, col] = ind2sub([17,8], q);
+        normalITD(row,col,p) = itd_y(itd == round(condition.normal.avg.ITDm(row,col, cF_ind(p) )  )  );
+        normalILD(row,col,p) = ild_y(round(ild,2) == round(condition.normal.avg.ILDm(row,col,cF_ind(p)  ),2));
+        ruffcutITD(row,col,p) = itd_y(itd == round(condition.ohne.avg.ITDm(row,col, cF_ind(p) )  )  );
+        ruffcutILD(row,col,p) = ild_y(round(ild,2) == round(condition.ohne.avg.ILDm(row,col,cF_ind(p)  ),2));
+
+    end
+    
+end
+normaltun = (normalITD.*normalILD)*.01;
+ruffcuttun = (ruffcutITD.*ruffcutILD)*.01;
+
+%% model neuron tuning across itds and freqs, normal HRTF
+
+normalmap = figure;
+set(gcf,'Position',[100 100 900 300]);
+
+sgtitle('Normal', FontSize=20)
+
+
+for plotID = 2:length(freq)
+    subplot(1,length(freq), plotID-1)
+    imagesc(azimuths(6:12), elevations(2:6), normaltun(6:12,2:6,plotID)'); axis xy;  caxis([0 1]); colormap jet;
+    title([num2str(freq(plotID)/1000), ' kHz'], FontSize=16);
+    xticks([azimuths(6:3:12)])
+    yticks([elevations(2:2:6)])
+
+    ax = gca;
+    ax.FontSize = 16;
+    if plotID == 2
+        ylabel('Elevation (°)', FontSize=16)
+        e = text(-200,90, 'b');
+        e.FontSize = 20;
+    end
+    if plotID == 2
+        xlabel('Azimuth (°)', FontSize=16);
+
+    end
+end
+
+subplot(1,length(freq), length(freq))
+cc = imagesc(azimuths, elevations, normaltun(:,:,plotID)'); axis xy;  caxis([0 1]); colormap jet;
+c = colorbar('manual', 'Position', [.82 0.25 0.02, 0.5]);
+c.FontSize = 16;
+d = text(0, -90,'Norm. Spike Count');
+d.FontSize = 16; d.Rotation = 90; 
+set(cc,'Visible', 'off')
+set(get(cc,'Children'),'Visible','off');
+
+
+axis off
+%saveas(normalmap, 'F:\Various Files\My Papers\Shadron and Pena 2022\eLife Figures\normalmap.png')
+
+
+%% model neuron tuning across itds and freqs, ruff-removed HRTF
+
+ruffcutmap = figure;
+
+set(gcf,'Position',[100 100 900 300]);
+
+sgtitle('Ruff-Removed', FontSize=20)
+
+
+for plotID = 2:length(freq)
+    subplot(1,length(freq), plotID-1)
+    imagesc(azimuths(6:12), elevations(2:6), ruffcuttun(6:12,2:6,plotID)'); axis xy;  caxis([0 1]); colormap jet;
+    title([num2str(freq(plotID)/1000), ' kHz'], FontSize=16);
+    xticks([azimuths(6:3:12)])
+    yticks([elevations(2:2:6)])
+
+    ax = gca;
+    ax.FontSize = 16;
+    if plotID == 2
+        ylabel('Elevation (°)', FontSize=16)
+        e = text(-200,90, 'c');
+        e.FontSize = 20;
+    end
+    if plotID == 2
+        xlabel('Azimuth (°)', FontSize=16);
+
+    end
+end
+
+subplot(1,length(freq), length(freq))
+cc = imagesc(azimuths, elevations, ruffcuttun(:,:,plotID)'); axis xy;  caxis([0 1]); colormap jet;
+c = colorbar('manual', 'Position', [.82 0.25 0.02, 0.5]);
+c.FontSize = 16;
+d = text(0, -90,'Norm. Spike Count');
+d.FontSize = 16; d.Rotation = 90; 
+set(cc,'Visible', 'off')
+set(get(cc,'Children'),'Visible','off');
+
+axis off
+%saveas(ruffcutmap, 'F:\Various Files\My Papers\Shadron and Pena 2022\eLife Figures\ruffcutmap.png')
+
+%% Other subplots for Figure 5
+
+plotID = 5;
+n = 29;
+
+ITDmap = figure;
+set(gcf,'Position',[100 100 400 300]);
+imagesc(azimuths(6:12), elevations(2:6), condition.normal.avg.ITDm((6:12),(2:6),n)'); axis xy;  caxis([-150 150]);
+title(['Frequency: ', num2str(cF(n)/1000), ' kHz'], FontSize=16); 
+c = colorbar('XTick',[-150 0 150]); c.FontSize = 12; c.Label.FontSize = 18; colormap jet;
+ax = gca; xticks(azimuths(6:12)); yticks(elevation(2:6)); ax.FontSize = 14; ylabel('Elevation (°)'); xlabel('Azimuth (°)')
+set(c.XLabel,{'String','Rotation','Position'},{'ITD (μs)',90,[2 0]})
+%saveas(ITDmap, 'F:\Various Files\My Papers\Shadron and Pena 2022\Matlab\ITDmap.png')
+
+ILDmap = figure;
+set(gcf,'Position',[100 100 400 300]);
+imagesc(azimuths(6:12), elevations(2:6), condition.normal.avg.ILDm((6:12),(2:6),n)'); axis xy;  caxis([-16 16]);
+title(['Frequency: ', num2str(cF(n)/1000), ' kHz'], FontSize=16); 
+c = colorbar('XTick',[-15 0 15]); c.FontSize = 14; c.Label.FontSize = 18; colormap jet;
+ax = gca; xticks(azimuths(6:12)); yticks(elevation(2:6)); ax.FontSize = 14; ylabel('Elevation (°)'); xlabel('Azimuth (°)');
+set(c.XLabel,{'String','Rotation','Position'},{'ILD (dB)',90,[2 0]})
+%saveas(ILDmap, 'F:\Various Files\My Papers\Shadron and Pena 2022\Matlab\ILDmap.png')
+
+ITDcurve = figure;
+set(gcf,'Position',[100 100 200 150]);
+plot(itd, itd_y, LineWidth=3, Color = 'black');
+xlim([-200,200]); ax = gca; ax.FontSize = 12; ylabel('Spike Count'); xlabel('ITD (μs)')
+%saveas(ITDcurve, 'F:\Various Files\My Papers\Shadron and Pena 2022\Matlab\ITDcurve.png')
+
+ILDcurve = figure;
+set(gcf,'Position',[100 100 200 150]);
+plot(ild, ild_y, LineWidth=3, Color = 'black');
+xlim([-16,16]); ax = gca; ax.FontSize = 12; ylabel('Spike Count'); xlabel('ILD (dB)');
+%saveas(ILDcurve, 'F:\Various Files\My Papers\Shadron and Pena 2022\Matlab\ILDcurve.png')
+
+sixkhz = figure;
+set(gcf,'Position',[100 100 400 300]);
+imagesc(azimuths(6:12), elevations(2:6), normaltun(6:12,2:6,plotID)'); axis xy; caxis([0 1]); colormap jet;
+title('Spatial Tuning', FontSize=16); ax = gca; ax.FontSize = 16;
+c = colorbar; c.Label.FontSize = 16; 
+ylabel('Elevation (°)'); xlabel('Azimuth (°)'); xticks(azimuths(6:12)); yticks(elevation(2:6));
+set(c.XLabel,{'String','Rotation','Position'},{'Norm. Spike Count',90,[2.5 .5]})
+%saveas(sixkhz, 'F:\Various Files\My Papers\Shadron and Pena 2022\Matlab\spatialtun.png')

circmean.m

@@ -0,0 +1,27 @@
+function cm=circmean(ph)
+% Synopsis
+% cm=circmean(ph)
+% Desription:
+% cm is the circular mean or mean direction of the phase data in 
+% ph (angular data) in radians
+% ignores NaN
+% for ma matrix, calculates the mean for every column
+[rows,cols]=size(ph);
+if(rows==1 & cols>1) 
+    ph=ph';
+    [rows,cols]=size(ph);
+end;    
+for(c=1:cols)
+    C=sum(cos(ph(:,c)), 'omitnan');
+    S=sum(sin(ph(:,c)), 'omitnan');
+    t=atan(S/C);
+    if (C<0)
+        cm(c)=t+pi;
+    else
+        if (S<0)
+            cm(c)=t+2*pi;
+        else
+            cm(c)=t;
+        end;
+    end;
+end;

circstd.m

@@ -0,0 +1,7 @@
+function cstd=circstd(ph)
+% Synopsis
+% cstd=circstd(ph)
+% Desription:
+% cstd is the circular standarddeviation of the phase data in 
+% ph (angular data in radians) in radians
+cstd=sqrt(-2*log(1-circvar(ph)));

circvar.m

@@ -0,0 +1,15 @@
+function cvar=circvar(ph)
+% Synopsis
+% cm=circvar(ph)
+% Desription:
+% cm is the circular variance of the phase data in 
+% ph (angular data) in radians
+C=sum(cos(ph), 'omitnan');
+S=sum(sin(ph), 'omitnan');
+R=sqrt(S^2+C^2);
+l=find(not(isnan(ph)));
+if(length(l)>0);
+   cvar=1-R/length(l);
+else
+   cvar=NaN;
+end;

f_gsafecmap.m

@@ -0,0 +1,42 @@
+%f_gsafecmap.m
+function gcmap = f_gsafecmap(nEl)
+% This generates a colormap (gcmap) that renders well in grayscale.
+% This is a refinement of Carey Rappaport's colormap CMRmap,
+%
+% http://www.mathworks.com/matlabcentral/fileexchange/2662-cmrmap-m
+%
+% C. Rappaport, “A color map for effective black-and-white rendering
+% of color-scale images�, IEEE Antennas Propagat. Mag., vol. 44, no.
+% 3, pp. 94–96, Jun. 2002.
+%
+%
+% Optional Input: nEl - number of elements in colormap
+ 
+if nargin<1
+  nEl = 64;
+end
+ 
+CMRmap = [ ...
+0.00,0.00,0.00; ...
+0.15,0.15,0.50; ...
+0.30,0.15,0.75; ...
+0.60,0.20,0.50; ...
+1.00,0.25,0.15; ...
+0.90,0.50,0.00; ...
+0.90,0.75,0.10; ...
+0.90,0.90,0.50; ...
+1.00,1.00,1.00];
+ 
+x9 = linspace(0,1,9);
+xcmap = linspace(0,1,nEl);
+ 
+% resample to nEl values using splines
+for i = 1:3
+  gcmap(:,i) = interp1(x9,CMRmap(:,i),xcmap,'spline');
+end
+ 
+% maintain bounds
+gcmap(gcmap<0)=0;
+gcmap(gcmap>1)=1;
+ 
+end