alan@461: # Carfac.py - Cochlear filter model based on Dick Lyons work.  This material taken from his Hearing book (to be published)
alan@470: # Author: Al Strelzoff: strelz@mit.edu
alan@461: 
alan@470: # The point of this file is to extract some of the formulas from the book and practice with them so as to better understand the filters.
alan@470: # The file has been written and tested for python 2.7
alan@461: 
alan@470: from numpy import cos, sin, pi, e, real,imag,arctan2
alan@470: 
alan@470: 
alan@470: from pylab import figure, plot,loglog, title, axis, show
alan@461: 
alan@461: fs = 22050.0            # sampling rate
alan@461: 
alan@461: 
alan@461: # given a frequency f, return the ERB
alan@461: def ERB_Hz(f):
alan@461: #    Ref: Glasberg and Moore: Hearing Research, 47 (1990), 103-138
alan@461:     return 24.7 * (1.0 + 4.37 * f / 1000.0)
alan@461: 
alan@461: 
alan@461: # ERB parameters
alan@461: ERB_Q = 1000.0/(24.7*4.37)  # 9.2645
alan@461: ERB_break_freq = 1000/4.37  # 228.833
alan@461: 
alan@461: ERB_per_step = 0.3333        
alan@461: 
alan@461: # set up channels
alan@461: 
alan@461: first_pole_theta = .78 * pi # We start at the top frequency.
alan@461: pole_Hz = first_pole_theta * fs / (2.0*pi)  # frequency of top pole
alan@461: min_pole_Hz = 40.0    # bottom frequency
alan@461: 
alan@461: # set up the pole frequencies according to the above parameters
alan@461: pole_freqs = []     # empty list of pole frequencies to fill, zeroth will be the top
alan@461: while pole_Hz >  min_pole_Hz:
alan@461:     pole_Hz = pole_Hz - ERB_per_step * ERB_Hz(pole_Hz)
alan@461:     pole_freqs.append(pole_Hz)
alan@461:     
alan@461: n_ch = len(pole_freqs)      # n_ch is the number of channels or frequency steps
alan@461: print('num channels',n_ch)
alan@461: 
alan@461: # Now we have n_ch, the number of channels, so can make the array of filters by instantiating the filter class (see below)
alan@461: 
alan@461: # before we make the filters, let's plot the position of the frequencies and the values of ERB at each.
alan@461: 
alan@461: fscale = []
alan@461: erbs = []
alan@461: 
alan@470: figure(1)
alan@461: for i in range(n_ch):
alan@461:     
alan@461:     f = pole_freqs[i]        # the frequencies from the list
alan@461:     ERB = ERB_Hz(f)       # the ERB value at each frequency
alan@461:     fscale.append(f)
alan@461:     erbs.append(ERB)
alan@461:     
alan@461: # plot a verticle hash at each frequency:
alan@461:     u = []
alan@461:     v = []
alan@461:     for j in range(5):
alan@461:         u.append(f)
alan@461:         v.append(10.0 + float(j))
alan@461: 
alan@461:     plot(u,v)  
alan@461:     
alan@461: loglog(fscale,erbs)
alan@461: 
alan@461: title('ERB scale')
alan@461: 
alan@461: 
alan@461: 
alan@461: # This filter class includes some methods useful only in design.  They will not be used in run time implementation.
alan@461: #  From figure 14.3 in Dick Lyon's book.
alan@470: #  When translating to C++, this class will become a struct, and all the methods will be moved outside.
alan@461: 
alan@461: 
alan@461: #########################################################The Carfac filter class#################################################################################
alan@461: 
alan@461: # fixed parameters
alan@461: min_zeta = 0.12
alan@461: 
alan@461: class carfac():
alan@461:     
alan@461:     
alan@461:     # instantiate the class (in C++, the constructor)
alan@461:     def __init__(self,f):
alan@461:         
alan@461:         self.frequency = f
alan@461:        
alan@461:         theta = 2.0 * pi * f/fs
alan@461:         r =  1.0 - sin(theta) * min_zeta
alan@461:         a = r * cos(theta)
alan@461:         c = r * sin(theta)
alan@470:         h = sin(theta)
alan@470:         g = (1.0 - 2.0 * a + r ** 2)/(1.0 - 2.0 * a + h * c + r ** 2)
alan@470:       
alan@470:        
alan@461:         
alan@470:         self.gh = g*h       # no need to repeat in real time
alan@461:         
alan@461:         # make all parameters properties of the class 
alan@461:         self.a = a
alan@461:         self.c = c
alan@461:         self.r = r
alan@461:         self.theta = theta
alan@461:         self.h = h
alan@461:         self.g = g
alan@461:         
alan@461:         
alan@461:         # the two storage elements.  Referring to diagram 14.3 on p.263, z2 is the upper storage register, z1, the lower
alan@461:         self.z1 = 0.0    
alan@461:         self.z2 = 0.0    
alan@461:       
alan@461:         
alan@461:         # frequency response of this filter
alan@461:         self.H = []
alan@461:         
alan@461:        
alan@461:         
alan@461:         # the total  frequency magnitude of this filter including all the filters in front of this one
alan@461:         self.HT = []       # this list will be filled by multiplying all the H's ahead of it together with its own (H)
alan@461:         
alan@461:         
alan@461:     
alan@461:     # execute one clock tick.  Take in one input and output one result. Execution semantics taken from fig. 14.3
alan@461:     #  This execution model is not tested in this file.  Here for reference.  See the file Exec.py for testing this execution model.  This is the main run time method.
alan@461:     def input(self,X):
alan@461: 
alan@461:         # recover the class definitions of these variables.  These statements below take up zero time at execution since they are just compiler declarations.
alan@470:         # computation below is organized as some loads, followed by 3 2x2 multiply accumulates
alan@470:         #  Note: this function is not exercised in this file and is here only for reference
alan@461:         
alan@461:         a = self.a
alan@461:         c = self.c
alan@461:         g = self.g
alan@470:         gh = self.gh
alan@461:         z1 = self.z1                    # z1 is the lower storage in fig. 14.3
alan@461:         z2 = self.z2
alan@461:         
alan@461:         # calculate what the next value of z1 will be, but don't overwrite current value yet.   
alan@470:         next_z1 = (a * z1) + (c * z2)   # Note:  it is a 2 element multiply accumulate. compute first so as not to have to do twice.
alan@461:         # the output Y
alan@470:         Y = g * X + gh * next_z1       # Note: organized as a 2 element multiply accumulate.
alan@461:         
alan@461:         #stores
alan@470:         self.z2 = X + (a * z2) - (c * z1)        #Note: this is a 2 element multiply accumulate
alan@470:         self.z1 = next_z1
alan@461:         
alan@461:         return Y                        # The output
alan@461:     
alan@461:     # complex frequency response of this filter at frequency w.  That is, what it contributes to the cascade
alan@461:     # this method is used for test only.  It finds the frequency magnitude.  Not included in run time filter class.
alan@461:     def Hw(self,w):
alan@461: 
alan@461:         a = self.a
alan@461:         c = self.c
alan@461:         g = self.g
alan@461:         h = self.h
alan@461:         r = self.r
alan@461:         z = e ** (complex(0,w))         # w is in radians so this is z = exp(jw)
alan@470:         return g * (1.0 + (h*c*z)/(z**2 - 2.0*a*z + r**2 ))      
alan@461:         
alan@461:     
alan@461:     
alan@461: ######################################################End of Carfac filter class########################################################################
alan@461:     
alan@461: # instantiate the filters
alan@461: 
alan@461: # n_ch is the number of filters as determined above
alan@461: 
alan@461: Filters = []    # the list of all filters, the zeroth is the top frequency
alan@461: for i in range(n_ch):
alan@461:     f = pole_freqs[i]
alan@461:     filter = carfac(f)  # note: get the correct parameters for r and h from Dick's matlab script.  Load them here from a table.
alan@461:     Filters.append(filter)
alan@461:     
alan@461: 
alan@461: 
alan@461: # sweep parameters
alan@461: steps = 1000
alan@461: 
alan@470: figure(2)   
alan@470: title('CarFac individual filter frequency response')
alan@470: # note:  the  scales are linear, not logrithmic as in the book
alan@461: 
alan@461: 
alan@461: for i in range(n_ch):
alan@461:     filter = Filters[i]
alan@461:     # plotting arrays
alan@461:     u = []
alan@461:     v = []
alan@461:     # calculate the frequency magnitude by stepping the frequency in radians
alan@461:     for j in range(steps):
alan@461:         
alan@470:         w = pi * float(j)/steps     
alan@461:         u.append(w)
alan@461:         mag = filter.Hw(w)      # freq mag at freq w
alan@461:         filter.H.append(mag)        # save for later use
alan@461:         filter.HT.append(mag)       # will be total response of cascade to this point after we do the multiplication in a step below
alan@470:         v.append(abs(mag))           # y plotting axis
alan@461:         
alan@461:       
alan@461:     plot(u,v)
alan@461: 
alan@461: 
alan@461: # calculate the phase response of the same group of  filters
alan@461: figure(3)
alan@470: title('Carfac individual filter Phase lag')
alan@461: 
alan@461: 
alan@461: for i in range(n_ch):
alan@461:     filter = Filters[i]
alan@461:     
alan@461:     u = []
alan@461:     v = []
alan@461:     for j in range(steps):
alan@470:         x = pi * float(j)/steps
alan@461:         
alan@461:         u.append(x)
alan@461:        
alan@461:         mag = filter.H[j]
alan@470:         phase = arctan2(-imag(mag),-real(mag))  - pi    # this formula used to avoid wrap around 
alan@461:        
alan@461:         v.append(phase)           # y plotting axis 
alan@461:      
alan@461:     plot(u,v)
alan@470:     axis([0.0,pi,-3.0,0.05])
alan@461:     
alan@461:     
alan@470: # calulate and plot cascaded frequency response    
alan@461: 
alan@461: figure(4)   
alan@470: title('CarFac cascaded filter frequency response')
alan@461: 
alan@461: 
alan@461: for i in range(n_ch-1):
alan@461:    
alan@461:     filter = Filters[i]
alan@461:     next = Filters[i+1]
alan@461:     
alan@461:     
alan@461:     u = []
alan@461:     v = []
alan@461:     for j in range(steps):
alan@470:         w = pi * float(j)/steps    
alan@470:         u.append(w)
alan@461:         mag = filter.HT[j] * next.HT[j]
alan@470:         next.HT[j] = mag
alan@470:         v.append(abs(mag))
alan@461: 
alan@461:       
alan@461:     plot(u,v)
alan@461:    
alan@461:   
alan@461: 
alan@461: # calculate and plot the phase responses of the cascaded filters
alan@461: 
alan@470: figure(5)
alan@470: title('Carfac cascaded filter Phase lag')
alan@461: 
alan@461: 
alan@461: for i in range(n_ch):
alan@470: 
alan@461:     filter = Filters[i]
alan@461:     
alan@461: 
alan@461:     u = []
alan@470:     v = []          # store list of phases
alan@470:     w = []          # second copy of phases needed for phase unwrapping
alan@461:     for j in range(steps):
alan@470:         x = pi * float(j)/steps
alan@461:         
alan@461:         u.append(x)
alan@461:         mag = filter.HT[j]
alan@470:         a = imag(mag)
alan@470:         b = real(mag)
alan@470:         phase = arctan2(a,b) 
alan@461:         
alan@470:         v.append(phase)  
alan@470:         w.append(phase)
alan@461:         
alan@470:     # unwrap the phase
alan@461: 
alan@461: 
alan@470:     for i in range(1,len(v)):
alan@470:         diff = v[i]-v[i-1]
alan@470:         if diff > pi:
alan@470:             for j in range(i,len(w)):
alan@470:                 w[j] -= 2.0 * pi
alan@470:         elif diff < -pi:
alan@470:             for j in range(i,len(w)):
alan@470:                 w[j] += 2.0 * pi
alan@470:                 
alan@470:         else: continue
alan@470:    
alan@470:      
alan@470:     plot(u,w)
alan@470:     axis([0.0,pi,-25.0,0.05])
alan@461: show()
alan@461: