Mercurial > hg > aimc

comparison trunk/matlab/bmm/carfac/Carfac.py @ 706:f8e90b5d85fd tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Delete CARFAC code from this repository. It has been moved to https://github.com/google/carfac Please email me with your github username to get access. I've also created a new mailing list to discuss CARFAC development: https://groups.google.com/forum/#!forum/carfac-dev
author ronw@google.com
date Thu, 18 Jul 2013 20:56:51 +0000
parents be2b68ced23d
children
comparison
equal deleted inserted replaced
705:be2b68ced23d 706:f8e90b5d85fd
1 # Carfac.py - Cochlear filter model based on Dick Lyons work. This material taken from his Hearing book (to be published)
2 # Author: Al Strelzoff: strelz@mit.edu
3
4 # 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.
5 # The file has been written and tested for python 2.7
6
7 from numpy import cos, sin, pi, e, real,imag,arctan2,log10
8
9
10 from pylab import figure, plot,loglog, title, axis, show
11
12 fs = 22050.0 # sampling rate
13
14
15 # given a frequency f, return the ERB
16 def ERB_Hz(f):
17 # Ref: Glasberg and Moore: Hearing Research, 47 (1990), 103-138
18 return 24.7 * (1.0 + 4.37 * f / 1000.0)
19
20
21 # ERB parameters
22 ERB_Q = 1000.0/(24.7*4.37) # 9.2645
23 ERB_break_freq = 1000/4.37 # 228.833
24
25 ERB_per_step = 0.3333
26
27 # set up channels
28
29 first_pole_theta = .78 * pi # We start at the top frequency.
30 pole_Hz = first_pole_theta * fs / (2.0*pi) # frequency of top pole
31 min_pole_Hz = 40.0 # bottom frequency
32
33 # set up the pole frequencies according to the above parameters
34 pole_freqs = [] # empty list of pole frequencies to fill, zeroth will be the top
35 while pole_Hz > min_pole_Hz:
36 pole_Hz = pole_Hz - ERB_per_step * ERB_Hz(pole_Hz)
37 pole_freqs.append(pole_Hz)
38
39 n_ch = len(pole_freqs) # n_ch is the number of channels or frequency steps
40 print('num channels',n_ch)
41
42 # Now we have n_ch, the number of channels, so can make the array of filters by instantiating the filter class (see below)
43
44 # before we make the filters, let's plot the position of the frequencies and the values of ERB at each.
45
46 fscale = []
47 erbs = []
48
49 figure(1)
50 for i in range(n_ch):
51
52 f = pole_freqs[i] # the frequencies from the list
53 ERB = ERB_Hz(f) # the ERB value at each frequency
54 fscale.append(f)
55 erbs.append(ERB)
56
57 # plot a verticle hash at each frequency:
58 u = []
59 v = []
60 for j in range(5):
61 u.append(f)
62 v.append(10.0 + float(j))
63
64 plot(u,v)
65
66 loglog(fscale,erbs)
67
68 title('ERB scale')
69
70
71
72 # This filter class includes some methods useful only in design. They will not be used in run time implementation.
73 # From figure 14.3 in Dick Lyon's book.
74 # When translating to C++, this class will become a struct, and all the methods will be moved outside.
75
76
77 #########################################################The Carfac filter class#################################################################################
78
79 # fixed parameters
80 min_zeta = 0.12
81
82 class carfac():
83
84
85 # instantiate the class (in C++, the constructor)
86 def __init__(self,f):
87
88 self.frequency = f
89
90 theta = 2.0 * pi * f/fs
91 r = 1.0 - sin(theta) * min_zeta
92 a = r * cos(theta)
93 c = r * sin(theta)
94 h = sin(theta)
95 g = (1.0 - 2.0 * a + r ** 2)/(1.0 - 2.0 * a + h * c + r ** 2)
96
97
98
99 self.gh = g*h # no need to repeat in real time
100
101 # make all parameters properties of the class
102 self.a = a
103 self.c = c
104 self.r = r
105 self.theta = theta
106 self.h = h
107 self.g = g
108
109
110 # the two storage elements. Referring to diagram 14.3 on p.263, z2 is the upper storage register, z1, the lower
111 self.z1 = 0.0
112 self.z2 = 0.0
113
114
115 # frequency response of this filter
116 self.H = []
117
118
119
120 # the total frequency magnitude of this filter including all the filters in front of this one
121 self.HT = [] # this list will be filled by multiplying all the H's ahead of it together with its own (H)
122
123
124
125 # execute one clock tick. Take in one input and output one result. Execution semantics taken from fig. 14.3
126 # 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.
127 def input(self,X):
128
129 # recover the class definitions of these variables. These statements below take up zero time at execution since they are just compiler declarations.
130 # computation below is organized as some loads, followed by 3 2x2 multiply accumulates
131 # Note: this function is not exercised in this file and is here only for reference
132
133 a = self.a
134 c = self.c
135 g = self.g
136 gh = self.gh
137 z1 = self.z1 # z1 is the lower storage in fig. 14.3
138 z2 = self.z2
139
140 # calculate what the next value of z1 will be, but don't overwrite current value yet.
141 next_z1 = (a * z1) + (c * z2) # Note: it is a 2 element multiply accumulate. compute first so as not to have to do twice.
142 # the output Y
143 Y = g * X + gh * next_z1 # Note: organized as a 2 element multiply accumulate.
144
145 #stores
146 self.z2 = X + (a * z2) - (c * z1) #Note: this is a 2 element multiply accumulate
147 self.z1 = next_z1
148
149 return Y # The output
150
151 # complex frequency response of this filter at frequency w. That is, what it contributes to the cascade
152 # this method is used for test only. It finds the frequency magnitude. Not included in run time filter class.
153 def Hw(self,w):
154
155 a = self.a
156 c = self.c
157 g = self.g
158 h = self.h
159 r = self.r
160 z = e ** (complex(0,w)) # w is in radians so this is z = exp(jw)
161 return g * (1.0 + (h*c*z)/(z**2 - 2.0*a*z + r**2 ))
162
163
164
165 ######################################################End of Carfac filter class########################################################################
166
167
168 # instantiate the filters
169
170 # n_ch is the number of filters as determined above
171
172 Filters = [] # the list of all filters, the zeroth is the top frequency
173 for i in range(n_ch):
174 f = pole_freqs[i]
175 filter = carfac(f) # note: get the correct parameters for r and h from Dick's matlab script. Load them here from a table.
176 Filters.append(filter)
177
178
179
180 # sweep parameters
181 steps = 1000
182
183 figure(2)
184 title('CarFac individual filter frequency response')
185 # note: the scales are linear, not logrithmic as in the book
186
187
188 for i in range(n_ch):
189 filter = Filters[i]
190 # plotting arrays
191 u = []
192 v = []
193 # calculate the frequency magnitude by stepping the frequency in radians
194 for j in range(steps):
195
196 w = pi * float(j)/steps
197 u.append(w)
198 mag = filter.Hw(w) # freq mag at freq w
199 filter.H.append(mag) # save for later use
200 filter.HT.append(mag) # will be total response of cascade to this point after we do the multiplication in a step below
201 v.append(abs(mag)) # y plotting axis
202
203
204 plot(u,v)
205
206
207 # calculate the phase response of the same group of filters
208 figure(3)
209 title('Carfac individual filter Phase lag')
210
211
212 for i in range(n_ch):
213 filter = Filters[i]
214
215 u = []
216 v = []
217 for j in range(steps):
218 x = pi * float(j)/steps
219
220 u.append(x)
221
222 mag = filter.H[j]
223 phase = arctan2(-imag(mag),-real(mag)) - pi # this formula used to avoid wrap around
224
225 v.append(phase) # y plotting axis
226
227 plot(u,v)
228 axis([0.0,pi,-3.0,0.05])
229
230
231 # calulate and plot cascaded frequency response
232
233 figure(4)
234 title('CarFac cascaded filter frequency response')
235
236
237 for i in range(n_ch-1):
238
239 filter = Filters[i]
240 next = Filters[i+1]
241
242
243 u = []
244 v = []
245 for j in range(steps):
246 w = pi * float(j)/steps
247 u.append(w)
248 mag = filter.HT[j] * next.HT[j]
249 next.HT[j] = mag
250 v.append(abs(mag))
251
252
253 plot(u,v)
254
255
256
257 # calculate and plot the phase responses of the cascaded filters
258
259 figure(5)
260 title('Carfac cascaded filter Phase lag')
261
262 for i in range(n_ch):
263
264 filter = Filters[i]
265
266
267 u = []
268 v = [] # store list of phases
269 w = [] # second copy of phases needed for phase unwrapping
270 for j in range(steps):
271 x = pi * float(j)/steps
272
273 u.append(x)
274 mag = filter.HT[j]
275 a = imag(mag)
276 b = real(mag)
277 phase = arctan2(a,b)
278
279 v.append(phase)
280 w.append(phase)
281
282 # unwrap the phase
283
284
285 for i in range(1,len(v)):
286 diff = v[i]-v[i-1]
287 if diff > pi:
288 for j in range(i,len(w)):
289 w[j] -= 2.0 * pi
290 elif diff < -pi:
291 for j in range(i,len(w)):
292 w[j] += 2.0 * pi
293
294 else: continue
295
296 # convert delay to cycles
297 for i in range(len(w)):
298 w[i] /= 2.0 * pi
299
300
301 plot(u,w)
302 axis([0.0,pi,-5.0,0.05])
303
304 # calculate group delay as cascaded lag/filter center frequency
305
306 figure(6)
307 title('Carfac group delay')
308
309 for i in range(n_ch):
310
311 filter = Filters[i]
312
313
314 u = []
315 v = [] # store list of phases
316 w = [] # second copy of phases needed for phase unwrapping
317 for j in range(1,steps):
318 x = pi * float(j)/steps
319
320 u.append(x)
321 mag = filter.HT[j]
322 a = imag(mag)
323 b = real(mag)
324 phase = arctan2(a,b)
325
326 v.append(phase)
327 w.append(phase)
328
329
330 for i in range(1,len(v)):
331 diff = v[i]-v[i-1]
332 if diff > pi:
333 for j in range(i,len(w)):
334 w[j] -= 2.0 * pi
335 elif diff < -pi:
336 for j in range(i,len(w)):
337 w[j] += 2.0 * pi
338
339 else: continue
340
341 # calculate the group delay from: GroupDelay(n)=-[((ph(n+1)-ph(n-1))/(w(n+1)-w(n-1))]/2pi # w in radians
342 # or gd[n] = - ((delta ph)/(delta w))/2pi
343 # delta ph = w[n+1] - w[n-1]
344 # delta w = pi* fs/steps # we divide up the frequency range fs into fractions of pi radians.
345 # gd[n] = - steps(w[n+1] - w[n-1])/fs
346
347 A = -(float(steps))/(2 *fs) # note. could have calculated this once outside the filter loop
348
349 # gd[n] = A *(w[n+1] - w[n-1])
350
351 gd = [] # in radians
352 gd.append(0.0) # pad to equalize size of u array, end points are then meaningless
353 for n in range(1,len(w)-1):
354 gd.append(A * (w[n+1] - w[n-1]))
355
356 gd.append(0.0)
357
358 plot(u,gd)
359 axis([0.1,pi,-0.01,0.05])
360
361 show()
362