|
c@0
|
1 """
|
|
c@0
|
2 A module used for auditory analysis.
|
|
c@0
|
3
|
|
c@0
|
4 Models currently implemented:
|
|
c@0
|
5 * Frequency-modulation analysis model based on the human auditory system.
|
|
c@0
|
6
|
|
c@0
|
7 Model implementations in progress:
|
|
c@0
|
8 * Glasberg and Moore model
|
|
c@0
|
9
|
|
c@0
|
10 Packaged dependencies:
|
|
c@0
|
11 * utils.py and/or utils.pyc
|
|
c@0
|
12 * erb.dat
|
|
c@0
|
13 * outMidFir.dat
|
|
c@0
|
14
|
|
c@0
|
15 External dependencies:
|
|
c@0
|
16 * scipy
|
|
c@0
|
17 * numpy
|
|
c@0
|
18 * copy
|
|
c@0
|
19 * matplotlib
|
|
c@0
|
20 """
|
|
c@0
|
21
|
|
c@0
|
22 import utils
|
|
c@0
|
23 import scipy.signal as sp
|
|
c@0
|
24 import numpy as np
|
|
c@0
|
25 from copy import deepcopy
|
|
c@0
|
26 import matplotlib.pyplot as plt
|
|
c@0
|
27
|
|
c@0
|
28 def get_partial_loudness():
|
|
c@0
|
29 """
|
|
c@0
|
30 A function to calculate the partial loudness of an excitation pattern at given ERB.
|
|
c@0
|
31
|
|
c@0
|
32 TO DO
|
|
c@0
|
33 """
|
|
c@0
|
34
|
|
c@0
|
35
|
|
c@0
|
36 return
|
|
c@0
|
37
|
|
c@0
|
38 def get_excitation_i(input, fs, SPL, rectify=False):
|
|
c@0
|
39 """
|
|
c@0
|
40 A function to calculate the excitation intensity of the input signal.
|
|
c@0
|
41
|
|
c@0
|
42 Parameters:
|
|
c@0
|
43 * input (type: array-like matrix of floats) - signal normalised to an amplitude range of -1 to 1. (Required)
|
|
c@0
|
44 * fs (type: numerical) - sample frequency of the signal, input. (Required)
|
|
c@0
|
45 * SPL (type: double) - the sound pressure level (SPL) at 0 dBFS, e.g., the SPL of a sine
|
|
c@0
|
46 wave with peaks at amplitude 1 and troughs at amplitude -1. (Required)
|
|
c@0
|
47 * rectify (type: boolean) - Specifies whether to include half wave rectification, modelling the direction
|
|
c@0
|
48 of that the cochlear nerves vibrate.
|
|
c@0
|
49 True to include, False to ignore. (Optional; Default = False)
|
|
c@0
|
50
|
|
c@0
|
51 Returns:
|
|
c@0
|
52 * gtfs (type: numpy array of floats) - array with size ((39,) + np.shape(input)) containing the excitation
|
|
c@0
|
53 pattern (in sound intensity) for each ERB of input signal. The excitation pattern for the nth ERB can
|
|
c@0
|
54 be accessed with gtfs[n].
|
|
c@0
|
55 """
|
|
c@0
|
56
|
|
c@0
|
57 input = np.array(input)
|
|
c@0
|
58 inputOMFir = outMidFir(input)
|
|
c@0
|
59 inputPa = v_Pascal(inputOMFir, SPL)
|
|
c@0
|
60 gtfs = gamma_tone_filter(inputPa, fs)
|
|
c@0
|
61 if (rectify):
|
|
c@0
|
62 gtfs = half_rectification(gtfs)
|
|
c@0
|
63 gtfs = pa_i(gtfs)
|
|
c@0
|
64 gtfs = at_normalise(gtfs)
|
|
c@0
|
65
|
|
c@0
|
66 return gtfs
|
|
c@0
|
67
|
|
c@0
|
68 def plot_excitation_tf(fs = 44100, outMidFilt = True, xscale = 'log', yscale = 'log'):
|
|
c@0
|
69 """
|
|
c@0
|
70 A function that plots the transfer function of the outer middle ear and each gammatone filter.
|
|
c@0
|
71
|
|
c@0
|
72 Parameters:
|
|
c@0
|
73 * fs (type: numerical) - the sampling frequency of the signal. (Optional; Default = 44100)
|
|
c@0
|
74 * outMidFilt (type: boolean) - filter the signal by the outer and middle ear FIR filter. (Optional; Default = True)
|
|
c@0
|
75 * xscale (type: string) - the scale of the frequency axis. Values are 'log' or 'linear'. (Optional; Default = 'log')
|
|
c@0
|
76 * yscale (type: string) - the scale of the amplitude axis. Values are 'log' or 'linear'. (Optional; Default = 'log')
|
|
c@0
|
77
|
|
c@0
|
78 Returns:
|
|
c@0
|
79 * y (type: numpy array of floats) - array with size ((39,np.ceil(fs))) containing the impulse response of the
|
|
c@0
|
80 signal at each gammatone filter (and optionally, including the outer middle ear response). The response at the
|
|
c@0
|
81 nth gammatone filter can be accessed by y[n].
|
|
c@0
|
82 """
|
|
c@0
|
83
|
|
c@0
|
84 input = np.zeros(np.ceil(fs))
|
|
c@0
|
85 input[0] = 1
|
|
c@0
|
86 if(outMidFilt): input = outMidFir(input)
|
|
c@0
|
87 y = gamma_tone_filter(input, fs)
|
|
c@0
|
88
|
|
c@0
|
89 for i in range(np.shape(y)[0]):
|
|
c@0
|
90 utils.plot_fft(y[i],xscale, yscale, False)
|
|
c@0
|
91
|
|
c@0
|
92 plt.show()
|
|
c@0
|
93
|
|
c@0
|
94 return y
|
|
c@0
|
95
|
|
c@0
|
96 def get_modulation_i(input, fs):
|
|
c@0
|
97 """
|
|
c@0
|
98 A function to calculate the modulation intensity of the input intensity signal. The function implements a
|
|
c@0
|
99 filter bank of bandpass filters with cut off frequencies ranging from 0.25 to 16 Hz.
|
|
c@0
|
100
|
|
c@0
|
101 Parameters:
|
|
c@0
|
102 * input (type: array-like matrix of floats) - the input intensity signal.
|
|
c@0
|
103 E.g., use get_excitation_i() to obtain excitation intensity and use as input.
|
|
c@0
|
104 * fs (type: numerical) - sampling frequency of input signal
|
|
c@0
|
105
|
|
c@0
|
106 Returns:
|
|
c@0
|
107 * y (type: numpy array of floats) - array with size ((10,) + np.shape(input)) containing the
|
|
c@0
|
108 modulation intensity of the signal at each modulation filter. The modulation intensity for the nth filter can
|
|
c@0
|
109 be accessed with y[n].
|
|
c@0
|
110 """
|
|
c@0
|
111
|
|
c@0
|
112
|
|
c@0
|
113 input = np.array(input)
|
|
c@0
|
114 b = fir_antialias(fs)
|
|
c@0
|
115 input_lp = sp.lfilter(b,(1),input_fr)
|
|
c@0
|
116 input_ds = downsample(input_lp, fs)
|
|
c@0
|
117 fc = np.array(utils.exp_sequence(-2,4,10))
|
|
c@0
|
118 bw = fc/2
|
|
c@0
|
119 y = decomposition(input_ds, fs, fc, bw)
|
|
c@0
|
120
|
|
c@0
|
121 return y
|
|
c@0
|
122
|
|
c@0
|
123 def outMidFir(input):
|
|
c@0
|
124 """
|
|
c@0
|
125 A function to filter the input signal with the response of the outer and middle ear.
|
|
c@0
|
126
|
|
c@0
|
127 Parameters:
|
|
c@0
|
128 * input (type: array-like matrix of floats) - signal normalised to an amplitude range of -1 to 1. (Required)
|
|
c@0
|
129
|
|
c@0
|
130 Returns:
|
|
c@0
|
131 * y (type: numpy array of floats) - array with dimensions equal to the input signal filtered by the response of
|
|
c@0
|
132 the outer and middle ear.
|
|
c@0
|
133 """
|
|
c@0
|
134
|
|
c@0
|
135 input = np.array(input)
|
|
c@0
|
136 b = utils.load_outMidFir_coeff()
|
|
c@0
|
137 y = sp.lfilter(b, (1), input)
|
|
c@0
|
138
|
|
c@0
|
139 return y
|
|
c@0
|
140
|
|
c@0
|
141 def gamma_tone_filter(input, fs):
|
|
c@0
|
142 """
|
|
c@0
|
143 A function to filter to decompose the input signal into 39 different gammatones filtered signals modelling the ERBs.
|
|
c@0
|
144
|
|
c@0
|
145 Parameters:
|
|
c@0
|
146 * input (type: array-like matrix of floats) - signal normalised to an amplitude range of -1 to 1. (Required)
|
|
c@0
|
147 * fs (type: numerical) - sample frequency of the signal, input. (Required)
|
|
c@0
|
148
|
|
c@0
|
149 Returns:
|
|
c@0
|
150 * y (type: numpy array of floats) - array with size ((39),np.shape(input)) containing the impulse response of the
|
|
c@0
|
151 signal at each gammatone filter. The response at the nth gammatone filter can be accessed by y[n].
|
|
c@0
|
152 """
|
|
c@0
|
153
|
|
c@0
|
154 input = np.array(input)
|
|
c@0
|
155 fc, bw = utils.load_erb_data()
|
|
c@0
|
156 y = decomposition(input,fs,fc,bw)
|
|
c@0
|
157
|
|
c@0
|
158 return y
|
|
c@0
|
159
|
|
c@0
|
160 def v_Pascal(input, SPL):
|
|
c@0
|
161 """
|
|
c@0
|
162 A function to convert a signal, normalised to an amplitude range of -1 to 1, to a signal represented in pressure (units: Pascal).
|
|
c@0
|
163
|
|
c@0
|
164 Parameters:
|
|
c@0
|
165 * input (type: array-like matrix of floats) - signal normalised to an amplitude range of -1 to 1. (Required)
|
|
c@0
|
166 * SPL (type: double) - the sound pressure level (SPL) at 0 dBFS, e.g., the SPL of a sine
|
|
c@0
|
167 wave with peaks at amplitude 1 and troughs at amplitude -1. (Required)
|
|
c@0
|
168
|
|
c@0
|
169 Returns:
|
|
c@0
|
170 * y (type: numpy array of floats) - array with dimensions equal to the input signal containing the input represented
|
|
c@0
|
171 as a pressure signal.
|
|
c@0
|
172 """
|
|
c@0
|
173
|
|
c@0
|
174 input = np.array(input)
|
|
c@0
|
175 y = np.sign(input)*(0.00002*10**(np.log10(np.abs(input))+dBFS/20))
|
|
c@0
|
176
|
|
c@0
|
177 return y
|
|
c@0
|
178
|
|
c@0
|
179 def pa_i(input, C=406):
|
|
c@0
|
180 """
|
|
c@0
|
181 A function to convert a pressure signal (unit: Pascal) to an intensity signal.
|
|
c@0
|
182
|
|
c@0
|
183 Parameters:
|
|
c@0
|
184 * input (type: array-like matrix of floats) - pressure signal (unit: Pascal) (Required)
|
|
c@0
|
185 * C (type: double) - the acoustic impedance of the air (Optional; Default = 406)
|
|
c@0
|
186
|
|
c@0
|
187 Returns:
|
|
c@0
|
188 * y (type: numpy array of floats) - array with dimensions equal to the input signal containing the input represented
|
|
c@0
|
189 as a pressure signal.
|
|
c@0
|
190 """
|
|
c@0
|
191
|
|
c@0
|
192 input = np.array(input)
|
|
c@0
|
193 y = (input**2) / C
|
|
c@0
|
194
|
|
c@0
|
195 return y
|
|
c@0
|
196
|
|
c@0
|
197 def at_normalise(input):
|
|
c@0
|
198 """
|
|
c@0
|
199 A function to normalise an intensity signal with the audibility threshold.
|
|
c@0
|
200
|
|
c@0
|
201 Parameters:
|
|
c@0
|
202 * input (type: array-like matrix of floats) - intensity signal (unit: Pascal) (Required)
|
|
c@0
|
203
|
|
c@0
|
204 Returns:
|
|
c@0
|
205 * y (type: numpy array of floats) - array with dimensions equal to the input signal containing the input normalised
|
|
c@0
|
206 with the audibility threshold.
|
|
c@0
|
207 """
|
|
c@0
|
208
|
|
c@0
|
209 input = np.array(input)
|
|
c@0
|
210 y = input / 1*(10**12)
|
|
c@0
|
211
|
|
c@0
|
212 return
|
|
c@0
|
213
|
|
c@0
|
214 def downsample(input, factor=100):
|
|
c@0
|
215 """
|
|
c@0
|
216 A function to downsample a signal, input, with sampling frequency, fs, by a downsample factor of factor.
|
|
c@0
|
217
|
|
c@0
|
218 NOTE: It is advised to use the fir_antialias() function before downsampling to remove any high frequencies
|
|
c@0
|
219 which would otherwise represented as low frequencies due to aliasing.
|
|
c@0
|
220
|
|
c@0
|
221 Parameters:
|
|
c@0
|
222 * input (type: array-like matrix of floats) - input signal. (Required)
|
|
c@0
|
223 * factor - downsample factor (Optional; Default = 100)
|
|
c@0
|
224
|
|
c@0
|
225 Returns:
|
|
c@0
|
226 * output (type: numpy array of floats) - array with outer dimensions equivalent to the to the input, and
|
|
c@0
|
227 inner dimension equal to np.floor(lenIn / factor) where lenIn is the length of the inner dimension.
|
|
c@0
|
228 """
|
|
c@0
|
229
|
|
c@0
|
230 input = np.array(input)
|
|
c@0
|
231 shapeIn = np.shape(input)
|
|
c@0
|
232 nDim = np.shape(shapeIn)
|
|
c@0
|
233 lenIn = shapeIn[nDim[0]-1]
|
|
c@0
|
234 lenOut = np.floor(lenIn / factor)
|
|
c@0
|
235 n = np.linspace(0,lenIn,lenOut, endpoint=False).astype(np.int)
|
|
c@0
|
236 output = input[...,n]
|
|
c@0
|
237
|
|
c@0
|
238 return output
|
|
c@0
|
239
|
|
c@0
|
240 def antialias_fir(fs, fc=100, order=64):
|
|
c@0
|
241 """
|
|
c@0
|
242 A function which returns the b coefficients for a lowpass fir filter with specified requirements.
|
|
c@0
|
243 Made specifically to remove aliasing when downsampling, but can be used for any application that
|
|
c@0
|
244 requires a lowpass filter.
|
|
c@0
|
245
|
|
c@0
|
246 Parameters:
|
|
c@0
|
247 * fs (type: numerical) - sampling frequency of the signal to be filtered (Required)
|
|
c@0
|
248 * fc (type: numerical) - unnormalised cut off frequency of the filter (Optional; Default = 100)
|
|
c@0
|
249 * order (type: numerical int) - the order of the filter (number of taps minus 1) (Optional; Default = 64)
|
|
c@0
|
250
|
|
c@0
|
251 Returns:
|
|
c@0
|
252 * b (type: numpy array of floats) - an array of size order + 1 (i.e. a coefficient for each tap)
|
|
c@0
|
253 """
|
|
c@0
|
254
|
|
c@0
|
255 nyquist = 0.5*fs
|
|
c@0
|
256 fcNorm = fc/nyquist
|
|
c@0
|
257 b = sp.firwin(order+1, fcNorm)
|
|
c@0
|
258
|
|
c@0
|
259 return b
|
|
c@0
|
260
|
|
c@0
|
261 def half_rectification(input):
|
|
c@0
|
262 """
|
|
c@0
|
263 A function which performs a half-wave rectification on the input signal.
|
|
c@0
|
264
|
|
c@0
|
265 Parameters:
|
|
c@0
|
266 * input (type: array-like matrix of floats) - input signal. (Required)
|
|
c@0
|
267
|
|
c@0
|
268 Returns:
|
|
c@0
|
269 * y (type: numpy array of floats) - an array with dimensions of input containing the half-wave rectification of
|
|
c@0
|
270 input.
|
|
c@0
|
271 """
|
|
c@0
|
272
|
|
c@0
|
273
|
|
c@0
|
274 y = np.array(input)
|
|
c@0
|
275 y[y<0] = 0
|
|
c@0
|
276
|
|
c@0
|
277 return y
|
|
c@0
|
278
|
|
c@0
|
279 def decomposition(input, fs, fc, bw, order=1024, verbose = False):
|
|
c@0
|
280 """
|
|
c@0
|
281 A function to run the input filter through a bandpass filter bank of length equal to the length of fc. Each
|
|
c@0
|
282 bandpass filter is designed by defining the centre frequencies, fc, and bandwidths, bw.
|
|
c@0
|
283
|
|
c@0
|
284 Parameters:
|
|
c@0
|
285 * input (type: array-like matrix of floats) - input signal. (Required)
|
|
c@0
|
286 * fs (type: numerical) - the sampling frequency of the input signal. (Required)
|
|
c@0
|
287 * fc (type: array-like vector of floats) - the centre off frequencies (unnormalised) of each bandpass filter.
|
|
c@0
|
288 The length of this vector determines the number of filters in the bank.
|
|
c@0
|
289 * bw (type: array-like vector of floats) - the bandwidths (unnormalised) of each bandpass filter. Must be equal
|
|
c@0
|
290 to or more than the length of fc. If the length is more, all elements exceeding the length of fc - 1 will be
|
|
c@0
|
291 ignored.
|
|
c@0
|
292 * order (type: numerical int) - the order of the filter (number of taps minus 1) (Optional; Default = 1024)
|
|
c@0
|
293 * verbose (type: boolean) - determines whether to display the current subroutine/progress of the procedure.
|
|
c@0
|
294 (Optional; Default = False)
|
|
c@0
|
295
|
|
c@0
|
296 Returns:
|
|
c@0
|
297 * y (type: numpy array of floats) - an array with inner dimensions equal to that of the input and outer dimension equal to
|
|
c@0
|
298 the length of fc (i.e. the number of bandpass filters in the bank) containing the outputs to each filter. The output
|
|
c@0
|
299 signal of the nth filter can be accessed using y[n].
|
|
c@0
|
300 """
|
|
c@0
|
301
|
|
c@0
|
302
|
|
c@0
|
303 input = np.array(input)
|
|
c@0
|
304 nFreqs = len(fc)
|
|
c@0
|
305 shape = (nFreqs,) + np.shape(input)
|
|
c@0
|
306 shape = shape[0:]
|
|
c@0
|
307 y = np.zeros(shape)
|
|
c@0
|
308
|
|
c@0
|
309 if(verbose): print "Running frequency decomposition."
|
|
c@0
|
310 for i in range(nFreqs):
|
|
c@0
|
311 if(verbose): print str(100.0*i/nFreqs) + "% complete."
|
|
c@0
|
312 low = fc[i]-bw[i]/2;
|
|
c@0
|
313 high = fc[i]+bw[i]/2;
|
|
c@0
|
314 if(verbose): print "Low: " + str(low) + "; High: " + str(high)
|
|
c@0
|
315 b = fir_bandpass(low, high, fs, order, verbose)
|
|
c@0
|
316 x = deepcopy(input)
|
|
c@0
|
317 y[i] = sp.lfilter(b,(1),x)
|
|
c@0
|
318
|
|
c@0
|
319 return y
|
|
c@0
|
320
|
|
c@0
|
321 def fir_bandpass(low, high, fs, order = 1024, verbose = False):
|
|
c@0
|
322 """
|
|
c@0
|
323 A function which returns the b coefficients for a bandpass fir filter with specified requirements.
|
|
c@0
|
324
|
|
c@0
|
325 Parameters:
|
|
c@0
|
326 * low - the lower cutoff frequency of the filter. (Required)
|
|
c@0
|
327 * high - the upper cutoff frequency of the filter. (Required)
|
|
c@0
|
328 * fs (type: numerical) - sampling frequency of the signal to be filtered. (Required)
|
|
c@0
|
329 * order (type: numerical int) - the order of the filter (number of taps minus 1) (Optional; Default = 1024)
|
|
c@0
|
330 * verbose (type: boolean) - determines whether to display the current progress (or info on the current subroutine)
|
|
c@0
|
331 of the procedure. (Optional; Default = False)
|
|
c@0
|
332
|
|
c@0
|
333 Returns:
|
|
c@0
|
334 * b (type: numpy array of floats) - an array of size order + 1 (i.e. a coefficient for each tap).
|
|
c@0
|
335 """
|
|
c@0
|
336
|
|
c@0
|
337 nyquist = 0.5*fs
|
|
c@0
|
338 lowNorm = low/nyquist
|
|
c@0
|
339 highNorm = high/nyquist
|
|
c@0
|
340 if(verbose): print "Low: " + str(lowNorm) + "; High: " + str(highNorm)
|
|
c@0
|
341 b = sp.firwin(order+1, [lowNorm, highNorm], pass_zero=False)
|
|
c@0
|
342
|
|
c@0
|
343 return b |
