Renamed files, added guinea pig and mouse

Author: Wayne Zhao

Cochlea.py

@@ -12,11 +12,13 @@ class Cochlea(FilterBank):
   '''
   Model of Cochlea (a subclass of FilterBank)
   '''
-  def __init__(self, types=None, species=None, CF0=20, length=20, l_factor=3.8, xs=None, cfs=None, rho=1.000, Ap=None, bp=None, Bu=None, gain_const=None, peak_magndb=None, Bpeak=None, fpeak=None, phiaccum=None, Nbeta=None, Nf=None, Qerb=None, ERBbeta=None, ERBf=None, Qn=None, Qn2=None, BWndBbeta=None, BWndBf=None, BWn2dBbeta=None, BWn2dBf=None, Sbeta=None, Sf=None, n=10, n2=3, betas=None, freqs=None):
+  def __init__(self, species=None, type=None, CF0=20, l_factor=3.8, length=20, xs=None, cfs=None, rho=1.000, Ap=None, bp=None, Bu=None, gain_const=None, peak_magndb=None, Bpeak=None, fpeak=None, phiaccum=None, Nbeta=None, Nf=None, Qerb=None, ERBbeta=None, ERBf=None, Qn=None, Qn2=None, BWndBbeta=None, BWndBf=None, BWn2dBbeta=None, BWn2dBf=None, Sbeta=None, Sf=None, n=10, n2=3, betas=None, freqs=None):
     '''
-    Initializes Cochlea. Most arguments are the same as for `FilterBank` object
+    Initializes Cochlea. Most arguments are the same as for `FilterBank` object.
 
     Attributes:
+      species: fast way to initialize Cochlea for certain species
+      num_filters: number of filters to model cochlea with. Default is 4.
       CF0: characteristic frequency (kHz) of base of cochlea
       length: length of cochlea (mm)
       l_factor: constant factor for cochlear model (mm)
@@ -28,30 +30,40 @@ def __init__(self, types=None, species=None, CF0=20, length=20, l_factor=3.8, xs
     '''
     if species is not None:
       CF0, l_factor, length = self._given_species(species)
-
     self.cochlea_length = length
+    self.cf = (lambda x: CF0*np.exp(-x/l_factor))
 
-    args = {'Ap':Ap, 'bp':bp, 'Bu':Bu, 'gain_const':gain_const, 'peak_magndb':peak_magndb,
-            'Bpeak':Bpeak, 'fpeak':fpeak, 'phiaccum':phiaccum, 'Nbeta':Nbeta, 'Nf':Nf, 'Qerb':Qerb, 'ERBbeta':ERBbeta, 'ERBf':ERBf, 'Qn':Qn, 'Qn2':Qn2, 'BWndBbeta':BWndBbeta, 'BWndBbeta':BWndBf, 'BWn2dBbeta':BWn2dBbeta, 'BWn2dBf':BWn2dBf, 'Sbeta':Sbeta, 'Sf':Sf, 'n':n, 'n2':n2, 'betas':betas, 'freqs':freqs}
-    num_filters = 1
-    for v in args.values():
-      if np.ndim(v) >= 1:
-        num_filters = max(num_filters, len(v))
+    type = 'P' if type is None else type
 
-    self.cf = (lambda x: CF0*np.exp(-x/l_factor))
-    if cfs is None:
-      if xs is None:
-        xs = np.linspace(0, length, num_filters)
-      cfs = [self.cf(x) for x in xs]
+    if species is not None:
+      xs = np.linspace(0, length, 4) # let user set?
+      cfs = self.cf(xs)
+      Ap = 0.3768 * np.exp(-0.1366 * cfs) # seems improbable
+      bp = [1., 1., 1., 1.]
+      Bu = 3.714 * np.exp(0.03123 * cfs) # same here
+      args = {'Ap':Ap, 'bp':bp, 'Bu':Bu, 'cf':cfs}
     else:
-      if xs is not None:
-        raise Exception('Please provide either only a list of all locations along cochlea or a list of all characteristic frequencies')
-      xs = [-np.log(c/CF0)*l_factor for c in cfs]
-    self.xs = xs
-    types = 'P' if types is None else types
+      args = {'Ap':Ap, 'bp':bp, 'Bu':Bu, 'gain_const':gain_const, 'peak_magndb':peak_magndb,
+              'Bpeak':Bpeak, 'fpeak':fpeak, 'phiaccum':phiaccum, 'Nbeta':Nbeta, 'Nf':Nf, 'Qerb':Qerb, 'ERBbeta':ERBbeta, 'ERBf':ERBf, 'Qn':Qn, 'Qn2':Qn2, 'BWndBbeta':BWndBbeta, 'BWndBbeta':BWndBf, 'BWn2dBbeta':BWn2dBbeta, 'BWn2dBf':BWn2dBf, 'Sbeta':Sbeta, 'Sf':Sf, 'n':n, 'n2':n2}
+      num_filters = 1
+      for v in args.values():
+        if np.ndim(v) >= 1:
+          num_filters = max(num_filters, len(v))
 
-    args['cf'] = cfs
-    super().__init__(topology='parallel', type=types, **args)
+      if cfs is None:
+        if xs is None:
+          xs = np.linspace(0, length, num_filters)
+        cfs = [self.cf(x) for x in xs]
+      else:
+        if xs is not None:
+          raise Exception('Please provide either only a list of all locations along cochlea or a list of all characteristic frequencies')
+        xs = [-np.log(cf/CF0)*l_factor for cf in cfs]
+
+      self.xs = xs
+      args['cf'] = cfs
+      args['betas'] = [betas for _ in range(num_filters)]
+      args['freqs'] = [freqs for _ in range(num_filters)]
+    super().__init__(topology='series', type=type, **args)
 
     apexmost_filter = self.filters[-1]
 
@@ -64,8 +76,6 @@ def __init__(self, types=None, species=None, CF0=20, length=20, l_factor=3.8, xs
     self.bp_fun = (lambda x: np.exp(np.interp(x, self.xs, np.log(np.array(bp)))))
     self.Bu_fun = (lambda x: np.exp(np.interp(x, self.xs, np.log(np.array(Bu)))))
 
-    # beta = w/2pi/CF(x)
-
     p = 1j*bp_apex - Ap_apex
     # k and Z both normalized to not depend on l
     self.wavenumber = (lambda beta: (beta/l_factor) * 2 * Bu_apex * (1j*beta + Ap_apex) / ((1j*beta - p)*(1j*beta - p.conjugate())))
@@ -76,18 +86,16 @@ def __init__(self, types=None, species=None, CF0=20, length=20, l_factor=3.8, xs
     self.Z = self.impedance
 
   @classmethod
-  def five_param(cls, types=None, aAp=None, bAp=None, bp=None, aBu=None, bBu=None, gain_const=None, peak_magndb=None, CF0=20, length=20, xs=None, rho=1.000):
+  def five_param(cls, type=None, aAp=None, bAp=None, bp=None, aBu=None, bBu=None, gain_const=None, peak_magndb=None, CF0=20, l_factor=3.8, length=20, xs=None, rho=1.000, betas=None, freqs=None):
     '''
     Five parameter parameterization of Cochlea from (Alkhairy 2019)
     '''
     if xs is None:
       xs = np.linspace(0, length, 4)
-    cf = (lambda x: CF0*np.exp(-x/3.8))
-    # print('cfx', cf(0))
-    # print(bAp, np.exp(bAp*cf(0)))
+    cf = (lambda x: CF0*np.exp(-x/l_factor))
     Ap_func = (lambda x: aAp*np.exp(bAp*cf(x)))
     Bu_func = (lambda x: aBu*np.exp(bBu*cf(x)))
-    cochlea = cls(types=types, Ap=[Ap_func(x) for x in xs], bp=bp, Bu=[Bu_func(x) for x in xs], gain_const=gain_const, peak_magndb=peak_magndb, CF0=CF0, length=length, xs=xs, rho=rho, species=None)
+    cochlea = cls(type=type, Ap=[Ap_func(x) for x in xs], bp=bp, Bu=[Bu_func(x) for x in xs], gain_const=gain_const, peak_magndb=peak_magndb, CF0=CF0, length=length, xs=xs, rho=rho, species=None, betas=betas, freqs=freqs)
     cochlea.Ap_fun = Ap_func
     cochlea.Bu_fun = Bu_func
     return cochlea
@@ -111,23 +119,24 @@ def filter_at_location(self, x_coord, gain_const=None, peak_magndb=None, type='P
     return Filter(type=type, Ap=self.Ap_fun(x_coord), bp=self.bp_fun(x_coord), Bu=self.Bu_fun(x_coord), gain_const=gain_const, peak_magndb=peak_magndb, cf=self.cf(x_coord))
 
   def _given_species(self, species):
-    if species is not None:
-      if species == 'chinchilla':
-        CF0 = 28.131
-        l_factor = 3.6649
-        length = 35
-      elif species == 'human':
-        CF0 = 20.823
-        l_factor = 7.2382
-        length = 20
-      # elif species == 'guinea pig' or species == 'guineapig':
-      # elif species == 'mouse'
-      xs = np.linspace(0, length, 5)[1:]
-      Ap = [0.3768 * np.exp(-0.1366 * CF0 * np.exp(-x/l_factor)) for x in xs] # this seems improbable
-      bp = [1., 1., 1., 1.]
-      Bu = [3.714 * np.exp(0.03123 * CF0 * np.exp(-x/l_factor)) for x in xs] # same here
+    if species == 'chinchilla': # Muller et al HR 2010
+      CF0 = 28.131
+      l_factor = 3.6649
+      length = 20
+    elif species == 'human':
+      CF0 = 20.823
+      l_factor = 7.2382
+      length = 35
+    elif species == 'guinea pig' or species == 'guineapig': # Tsuji and Liberman 1997 J. Comp. Neurol. 381:188-202
+      CF0 = 54.732
+      l_factor = 3.3178
+      length = 20
+    elif species == 'mouse':
+      CF0 = 71.130
+      l_factor = 1.8566
+      length = 5.13
     else:
-      l_factor = 3.8
+      raise Exception(f'"{species}" is an unsupported species')
 
     return (CF0, l_factor, length)
 
@@ -191,7 +200,7 @@ def plot_impedance(self, betas=None, setting='realimag', custom_title='Normalize
       show: `True` if plot is to be shown, `False` otherwise. Default is `True`.
       phase_in_rad: Show phase in radians if True or in cycles otherwise
     '''
-    # plot Z and normalized Z
+    # plot Z and normalized Z?
     if betas is None:
       betas = np.linspace(0.01, self.bp_apex*1.5, 10000)
     Zdata = self.Z_norm(betas)
@@ -235,7 +244,6 @@ def signal_response_heatmap(self, signal, len_xs=20, custom_title='Signal Heatma
     '''
     sigs = []
     cfs = []
-    # signal.plot()
     for x in np.linspace(0.025, self.cochlea_length, len_xs):
       fil = self.filter_at_location(x)
       cfs += [round(self.cf(x), 2)]

Filter.py

@@ -67,7 +67,7 @@ def __init__(self, ir=None, tf=None, coeffs=None, roots=None, type=None, Ap=None
     has_roots = (roots is not None)
     has_params = any(param is not None for param in [Ap, bp, Bu])
     has_chars = any(characteristic is not None for characteristic in [Bpeak, fpeak, phiaccum, Nbeta, Nf, Qerb, ERBbeta, ERBf, Qn, Qn2, BWndBbeta, BWndBf, BWn2dBbeta, BWn2dBf, Sbeta, Sf])
-    if sum([has_coeffs, has_roots, has_tf, has_params, has_chars]) != 1:
+    if sum([has_coeffs, has_roots, has_tf, has_ir, has_params, has_chars]) != 1:
       raise Exception('Exactly one filter representation should be used')
 
     self.in_terms_of_normalized = True
@@ -86,7 +86,13 @@ def __init__(self, ir=None, tf=None, coeffs=None, roots=None, type=None, Ap=None
     if Sf is not None: Sbeta = Sf * self.cf**2
     if Nf is not None: Nbeta = Nf * self.cf
 
-    if betas is None: betas = np.geomspace(0.01, 10, 10000) # is there a more adaptive way to pick betas if it is not provided
+    if betas is None:
+      maxbeta = 10
+      if bp is not None:
+        maxbeta = 3*bp
+      if Bpeak is not None:
+        maxbeta = 3*Bpeak
+      betas = np.linspace(0.01, maxbeta, 10000)
 
     # eventually refactor by making the __init__ of the three classes deal with the logic (which it actually already mostly does)
     if has_tf:
@@ -510,7 +516,7 @@ def pole_zero_plot(self, custom_title=None, show=True):
       ax.axhline(y=0, color='k', ls=':')
       ax.axvline(x=0, color='k', ls=':')
       ax.scatter([z.real for z in zeros], [z.imag for z in zeros], marker='o', facecolors='none', edgecolors='tab:orange')
-      ax.scatter([z.real for z in poles], [z.imag for z in poles], marker='x', edgecolors='tab:blue')
+      ax.scatter([z.real for z in poles], [z.imag for z in poles], marker='x', facecolors='tab:blue')
       ax.set_xlabel('Re(z)')
       ax.set_ylabel('Im(z)')
       plt.axis('equal')

FilterType.py

@@ -78,7 +78,19 @@ def __init__(self, uid, tf=None, ir=None, betas=None):
     if betas is None:
       betas = np.geomspace(0.01, 10, 10000)
     self.init_with_tf = (tf is not None)
-    chars = helpers.computedfiltercharacteristics(tfunc=tf, betas=betas)
+    if self.init_with_tf:
+      chars = helpers.computedfiltercharacteristics(tfunc=tf, betas=betas)
+    else:
+      # much faster than using approx_tf, even if it's not as accurate
+      siglen = len(betas)*2
+      samprate = betas[-1]*2
+      timestamps = np.arange(siglen)/samprate
+      irarr = [ir(t) for t in timestamps]
+      tfapprox = scipy.fft.rfft(irarr)
+      origbetas = scipy.fft.rfftfreq(siglen, d=1/samprate)
+      chars = helpers.computedfiltercharacteristics(tfunc=(lambda f: np.interp(f, origbetas, tfapprox)), betas=betas)
+
+      # np.interp
 
     if tf is not None:
       if ir is not None:

Signal.py

@@ -1,6 +1,8 @@
 import numpy as np
-import scipy as sp
 import scipy.fft
+import scipy.signal
+import scipy.integrate
+import scipy.io
 import matplotlib.pyplot as plt
 
 def tolerate(x, eps=1e-10):
@@ -48,15 +50,15 @@ def __init__(self, mode='t', data=[0 for _ in range(9)], fs=1, evenlen=True):
     # having self.length as length of self.mode_t keeps a bit more info (literally)
     if self.mode in ['t', 'ttilde']:
       self.mode_t = data
-      self.mode_f = sp.fft.rfft(data)
+      self.mode_f = scipy.fft.rfft(data)
       self.length = len(data)
     else:
       if self.mode == 'w':
         self.mode_f = data / 2 / np.pi
       else:
         self.mode_f = data
       self.length = 2*len(data)-1-int(evenlen) # -2 if evenlen is True and -1 if evenlen is False
-      self.mode_t = sp.fft.irfft(self.mode_f, self.length)
+      self.mode_t = scipy.fft.irfft(self.mode_f, self.length)
 
     self.func = (lambda t: self.at_time(t, tolerance=0)) # automatically in time domain, from_function for other domains
     self.timestamps = np.arange(self.length)/fs
@@ -67,7 +69,7 @@ def __init__(self, mode='t', data=[0 for _ in range(9)], fs=1, evenlen=True):
     self.mean = np.mean(self.mode_t)
     self.rms = np.mean([x**2 for x in self.mode_t])**0.5
 
-    self.analytic = sp.signal.hilbert(self.mode_t)
+    self.analytic = scipy.signal.hilbert(self.mode_t)
     self.hilbert = np.real(self.analytic)
     self.inst_phase = np.unwrap(np.angle(self.analytic))
 
@@ -91,7 +93,7 @@ def from_function(cls, mode='t', func=(lambda x: 0), fs=1, num_samples=9, evenle
     if mode in ['t', 'ttilde']:
       sample_points = np.arange(num_samples)/fs
     else:
-      sample_points = sp.fft.rfftfreq(2*num_samples-1-int(evenlen), 1/fs)
+      sample_points = scipy.fft.rfftfreq(2*num_samples-1-int(evenlen), 1/fs)
     S = cls(mode=mode, data=func(sample_points), fs=fs, evenlen=evenlen)
     S.func = func
     return S
@@ -125,17 +127,17 @@ def linear_chirp(cls, f_init=1, w_init=None, f_final=10, w_final=None, fs=1, num
     return cls.from_function(mode='t', func=(lambda t: np.cos(np.pi*t*(fi*(2-t/endtime) + ff*(t/endtime)))), fs=fs, num_samples=num_samples, evenlen=(num_samples%2==0))
     # evenlen is actually never used in this case since the mode is always 't', but just for completeness
 
-  # @classmethod
-  # def from_instantaneous_frequency(cls, freq_func=(lambda x: 0), freqs=None, init_phase=0, fs=1, num_samples=9):
-  #   if freqs is None:
-  #     # gaussian quadrature?
-  #     ws = [2*np.pi*freq_func(i/fs) for i in range(num_samples)]
-  #   else:
-  #     ws = [2*np.pi*f for f in freqs]
+  @classmethod
+  def from_instantaneous_frequency(cls, freq_func=(lambda x: 0), freqs=None, init_phase=0, fs=1, num_samples=9):
+    if freqs is None:
+      # gaussian quadrature?
+      ws = [2*np.pi*freq_func(i/fs) for i in range(num_samples)]
+    else:
+      ws = [2*np.pi*f for f in freqs]
 
-  #   phases = sp.integrate.cumulative_trapezoid(ws, dx=1/fs, initial=0)+init_phase
+    phases = scipy.integrate.cumulative_trapezoid(ws, dx=1/fs, initial=0)+init_phase
 
-  #   return cls(mode='t', data=np.cos(phases).tolist(), fs=fs)
+    return cls(mode='t', data=np.cos(phases).tolist(), fs=fs)
 
   def __iter__(self):
     self.__idx = -1
@@ -253,13 +255,15 @@ def at_time(self, t, tolerance=1e-10):
     if tolerate(num, eps=tolerance):
       return self.mode_t[round(num)%self.length]
 
-    f = [k/self.length for k in self.mode_f]
+    f = self.mode_f
+    freqstamps = self.freqstamps
 
-    tot = f[0].real
+    tot = np.real(f[0])
     for idx in range(1, len(f)):
-      tot += 2 * (f[idx] * np.exp(2j*np.pi*idx*t/self.length)).real
+      tot += 2 * np.real(f[idx] * np.exp(2j*np.pi*t*freqstamps[idx]))
     if self.length%2 == 0:
-      tot -= (f[idx] * np.exp(2j*np.pi*idx*t/self.length)).real
+      tot -= np.real(f[idx] * np.exp(2j*np.pi*t*freqstamps[idx]))
+    tot /= self.length
     if tolerate(tot, eps=tolerance):
       tot = round(tot)
     return tot
@@ -358,12 +362,12 @@ def moving_spectral_entropy(self, window_len=9):
     Hs = []
     for i in range(self.length-window_len+1):
       data = self.mode_t[i:i+9]
-      spectrum = abs(sp.fft.fft(data))**2
+      spectrum = abs(scipy.fft.fft(data))**2
       distribution = spectrum/sum(spectrum)
       Hs += [-sum(p*np.log(p) for p in distribution)/np.log(window_len)]
     return Hs
 
-  def spectrogram(self, win=sp.signal.windows.gaussian(30, std=5, sym=True), hop=1, mfft=200, custom_title='Spectrogram', show=True):
+  def spectrogram(self, win=scipy.signal.windows.gaussian(30, std=5, sym=True), hop=1, mfft=200, custom_title='Spectrogram', show=True):
     '''
     Generates spectrogram of Signal. Returns [SFFT data, bounds]. \
       Since the window has a small width, the resulting SFFT is \
@@ -377,7 +381,7 @@ def spectrogram(self, win=sp.signal.windows.gaussian(30, std=5, sym=True), hop=1
       show: `True` if plot is to be shown, `False` otherwise. Default is `True`.
     '''
     N = self.length
-    ft = sp.signal.ShortTimeFFT(np.array(win), hop, self.fs, mfft=mfft)
+    ft = scipy.signal.ShortTimeFFT(np.array(win), hop, self.fs, mfft=mfft)
     S = ft.spectrogram(np.array(self.mode_t))
     windowed_S = S[:, ft.lower_border_end[1]:ft.upper_border_begin(N)[1]]
     bounds = (0, ft.delta_t*len(windowed_S[0]), *ft.extent(N)[2:])
@@ -471,4 +475,4 @@ def as_sound(self, filename):
     Attributes:
       filename: Filename to save sound file under.
     '''
-    sp.io.wavfile.write(filename, self.fs, np.array(self.mode_t))
+    scipy.io.wavfile.write(filename, self.fs, np.array(self.mode_t))