1
2
3 import math
4 import numpy
5 from gmisclib import ortho_poly
6 from gmisclib import die
7 from gmisclib import gpkmisc
8
9 from gpk_voicing import gammatone as GT
10
11
13 """An instance of this class produces a window with a specified number of nodes.
14 So, you construct the class (telling it the window size), and that
15 gives you a class instance. Then, you use the class instance (by
16 calling the self.P method) it creates a window with the specified
17 number of nodes.
18 """
20 """@type n: int or None
21 @param n: how many points in the window
22 """
23 ortho_poly.ortho.__init__(self, n, x)
24
26 return numpy.square(1.0-numpy.square(self.x))
27
29 """@param nnodes: Now many zeros does the window have?
30 """
31 if nnodes%2 == 0:
32 wiggle = numpy.cos((math.pi*(nnodes//2))*self.x)
33 else:
34 wiggle = numpy.sin((math.pi*((nnodes+1)//2))*self.x)
35 return (1.0 + numpy.cos(math.pi*self.x)) * wiggle
36
37
38 _wincache = {}
39 _WIN_CACHE_SIZE = 20
40
42 """Make a window of length n with k nodes.
43 The window is normalized so that numpy.sum(w**norm)==1."""
44
45 assert n > 0
46 assert k >= 0
47 cachekey = (n,k,norm)
48 if cachekey in _wincache:
49 return _wincache[ cachekey ]
50
51 opn = cont_kernel_c(n)
52
53 w = opn.P(k)
54 if norm > 0:
55 numpy.divide(w, (numpy.absolute(w)**norm).sum()**(1.0/norm), w)
56 while len(_wincache) > _WIN_CACHE_SIZE:
57 _wincache.popitem()
58 _wincache[ cachekey ] = w
59 return w
60
61
62
71
72
74 """Make a window with width n,
75 except that it changes continuously when n
76 is a floating point number. The resulting window is
77 always an odd length.
78 """
79 assert n >= k
80 n2 = n/2.0
81 nu = 2*int(math.ceil(n2)) + 1
82 assert nu > 1
83 nl = 2*int(math.floor(n2)) + 1
84 assert nu==nl+2 or nu==nl
85 if nu==nl or nu<=k:
86 w = kernel(nu, k, norm=norm)
87 else:
88 f = n2 - math.floor(n2)
89 w = kernel(nu, k, norm=norm)
90 numpy.multiply(w, f, w)
91 wl = window(nl, k, norm=norm)
92 wul = w.shape[0]
93 dl = (wul-wl.shape[0])//2
94 assert dl>=0 and wul==wl.shape[0]+2*dl, "dl=%d wul=%d wll=%d" % (dl, wul, wl.shape[0])
95 numpy.multiply(wl, 1-f, wl)
96 numpy.add(w[dl:wul-dl], wl, w[dl:wul-dl])
97 if norm > 0:
98 numpy.divide(w, (numpy.absolute(w)**norm).sum()**(1.0/norm), w)
99 return w
100
101
103 """exp(a), except protected from underflows."""
104 MIN = -100
105 ok = numpy.greater(a, MIN)
106 return ok * numpy.exp(numpy.maximum(a, MIN))
107
108
109
110
112 """Find start and end points for a window.
113 @param t: desired window center (seconds).
114 @type t: float
115 @param window: window width (in samples).
116 @type window: int
117 @param dt: sampling interval (seconds).
118 @type dt: float
119 @param len_d: length of data array (samples).
120 @type len_d: int
121 @return (starting sample, end) where end=sample beyond the last
122 @rtype (int, int)
123 """
124 assert type(window) == type(1)
125 s = int(round(t/dt - 0.5*window))
126 if s < 0:
127 s = 0
128 elif s+window > len_d:
129 s = len_d - window
130 e = s + window
131 return (s, e)
132
133
135 np = n + max(10.0/cutoff_freq, n*0.01)
136 nx = near_win_size(min=np+1, tol=1+n//11, real=True)
137 return nx - nx%2
138
139
140
141 _fftworkcache = {}
143 """This estimates the work required to compute a FFT of size N,
144 using the Cooley-Tukey algorithm.
145 See http://en.wikipedia.org/wiki/Cooley-Tukey_FFT_algorithm
146 """
147 try:
148 return _fftworkcache[n]
149 except KeyError:
150 pass
151 f1 = gpkmisc.a_factor(n)
152 if f1 == n:
153 tmp = n**2
154 else:
155 nof1 = n/f1
156 tmp = f1*_fftwork(nof1) + n + nof1*_fftwork(f1)
157 _fftworkcache[n] = tmp
158 return tmp
159
160
162 """This is the work for a forward and reverse transform of
163 real data back to real data, using Num.FFT.rfft and Num.FFT.irfft
164 """
165 return _fftwork(n) + _fftwork(n-n%2)
166
167
168 -def near_win_size(near=None, tol=None, min=None, max=None, real=False):
169 """Find a window size near the specified size where
170 the FFT algorithm is fastest."""
171
172 if near is not None and tol is not None:
173 nmax = int(round(near + tol))
174 nmin = int(round(near - tol))
175 elif min is not None and max is not None:
176 nmax = int(round(max))
177 nmin = int(round(min))
178 elif min is not None and tol is not None:
179 nmin = int(round(min))
180 nmax = int(round(min + tol))
181 elif max is not None and tol is not None:
182 nmax = int(round(max))
183 nmin = int(round(max - tol))
184 else:
185 raise ValueError, "Needs near&tol, min&max, min&tol, or max&tol."
186
187 if real:
188 fftwork = _fftwork_real
189 else:
190 fftwork = _fftwork
191
192 besti = None
193 bestwork = None
194 Since_lim = 100
195 since = 0
196 for i in frommid_generator(nmin, nmax):
197 t = fftwork(i)
198 if bestwork is None or t < bestwork:
199 bestwork = t
200 besti = i
201 since = 0
202 else:
203
204
205
206
207 since += 1
208 if since > Since_lim:
209 break
210 return besti
211
212
214 """Generate numbers in [s,e] starting at the midpoint."""
215 mid = (s+e)//2
216 yield mid
217 for i in range(1, min(mid-s, e-mid)+1):
218 if mid+i <= e:
219 yield mid + i
220 if mid-i >= s:
221 yield mid - i
222
225
226
235
236
238 N = 100
239 w0 = window(N, 0)
240 w1 = window(N, 1)
241 w2 = window(N, 2)
242 w3 = window(N, 3)
243 w4 = window(N, 4)
244 w5 = window(N, 5)
245 for i in range(N):
246 print i, w0[i], w1[i], w2[i], w3[i], w4[i], w5[i]
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283 hipass_first_order = GT.hipass_first_order
284 lopass_first_order = GT.lopass_first_order
285
286
288 """A combinatio of a first-order high-pass filter
289 and a low-pass filter.
290 Cutoff freq is measured in cycles per point."""
291
292 assert cutoff_lp >= 0.0
293 assert cutoff_hp >= 0.0
294 assert len(d.shape) == 1
295 n = d.shape[0]
296 nx = _filter_win_size(n, 1.0/math.hypot(1.0/cutoff_lp, 1.0/cutoff_hp))
297
298 dx = numpy.fft.rfft(d, nx)
299 f = real_fft_freq(dx.shape[0])
300
301 frlp = 1j* (f / cutoff_lp)
302 frhp = 1j* (f / cutoff_hp)
303 ff = frhp / ((1 + frhp)*(1+frlp))
304 ff[0] = ff[0].real
305 ff[-1] = ff[-1].real
306 return numpy.fft.irfft(dx*ff, nx)[:n]
307
308
310 q = numpy.cos(numpy.arange(10000)*0.003*2*math.pi)
311 sq = numpy.sum(q*q)
312 q1 = hipass_first_order(q, 0.0003)
313 sq1 = numpy.sum(q1*q1)
314 assert 0.98*sq < sq1 < 0.995*sq
315 q4 = hipass_first_order(q, 0.003)
316 sq4 = numpy.sum(q4*q4)
317 assert 0.4*sq < sq4 < 0.6*sq
318 q5 = hipass_first_order(q, 0.03)
319 sq5 = numpy.sum(q5*q5)
320
321 assert 0.003 < sq5 < 0.03*sq
322
323
324
325
327 """Returns an array of indices whose frequencies are between f0 and f1
328 after you've called fft().
329 This will produce negative values where appropriate.
330 """
331 if hasattr(numpy.fft, 'fftfreq'):
332 return numpy.fft.fftfreq(n, d)
333 fidx = numpy.arange(n)
334 fneg = fidx - n
335 useneg = numpy.less(-fneg, fidx)
336 df = 0.5/(n*d)
337 return df * numpy.where(useneg, fneg, fidx)
338
339
341 """Returns an array of indices whose frequencies are between f0 and f1
342 after you've called real_fft().
343 """
344 raise RuntimeError, "Sorry! Not implemented"
345
346
348 """Returns the frequencies associated with an index, after you've
349 called real_fft().
350 @param d: The time interval between samples.
351 @type d: float
352 @param n: the size of the transformed array (i.e. frequency domain, complex)
353 @type n: int
354 @rtype: numpy.ndarray
355 @return: the frequency associated with each element in a Fourier Transform
356 """
357 fidx = numpy.arange(n)
358 df = 0.5/(n*d)
359 return df * fidx
360
361
363 """Returns the indices whose frequencies are between f0 and f1
364 after you've called real_fft().
365 """
366 df = 0.5/(n*d)
367 f0 = max(0.0, f0)
368 imin = int(math.ceil(f0/df))
369 imax = int(math.floor(f1/df))
370 return (imin, imax+1)
371
372
373
376
379
380
382 """Pad the data, subtract the average, and round
383 the edges of the window.
384 """
385 assert len(data.shape) == 1
386 m = data.shape[0]
387 if pad is None:
388 pad = 2*edgewidth + m//10
389 min_win_size = m + pad
390 n = near_win_size(min=min_win_size, tol=m//11 + 1, real=isreal)
391 d = numpy.zeros((n,))
392 avg = numpy.average(data)
393 d[:m] = data - avg
394 nedge = int(round(edgewidth/dt))
395 assert nedge <= m
396 ew = 0.5*(1.0-numpy.cos((math.pi/nedge)*(numpy.arange(nedge)+0.5)))
397 numpy.multiply(d[:nedge], ew, d[:nedge])
398 numpy.multiply(d[m-nedge:m], ew[::-1], d[m-nedge:m])
399 return d
400
401
403 """Pad the data, subtract the average, and round
404 the edges of the window.
405 """
406 assert len(data.shape) == 2
407 m = data.shape[0]
408 if pad is None:
409 pad = 2*edgewidth + m//10
410 min_win_size = m + pad
411 n = near_win_size(min=min_win_size, tol=m//11 + 1, real=isreal)
412 d = numpy.zeros((n,data.shape[1]))
413 avg = numpy.average(data, axis=0)
414 d[:m,:] = data - avg
415 nedge = int(round(edgewidth/dt))
416 assert nedge <= m
417 ew = 0.5*(1.0-numpy.cos((math.pi/nedge)*(numpy.arange(nedge)+0.5)))
418 numpy.multiply(d[:nedge], ew[:,numpy.newaxis], d[:nedge])
419 numpy.multiply(d[m-nedge:m], ew[::-1,numpy.newaxis], d[m-nedge:m])
420 return d
421
422
424 """A time-symmetric filter with a magnitude response that
425 is the same as the Butterworth filter.
426 Cutoff freq is measured in cycles per point."""
427 assert len(d.shape) == 1
428 if cutoff_freq<=0 or cutoff_freq>10.0:
429 die.warn('Silly value of cutoff_freq: %g cycles per point.' % cutoff_freq)
430 n = d.shape[0]
431 nx = _filter_win_size(n, cutoff_freq)
432 dx = numpy.fft.rfft(d, nx)
433 f = real_fft_freq(dx.shape[0])
434 ff = numpy.sqrt(1.0/(1 + (f/cutoff_freq)**(2*order)))
435 return numpy.fft.irfft(dx*ff, nx)[:n]
436
437
438
439
441 """A time-symmetric filter with a Gaussian impulse response.
442 Cutoff freq is measured in cycles per point.
443 """
444 assert len(d.shape) == 1
445 if cutoff_freq<=0 or cutoff_freq>10.0:
446 die.warn('Silly value of cutoff_freq: %g cycles per point.' % cutoff_freq)
447 n = d.shape[0]
448 nx = _filter_win_size(n, cutoff_freq)
449 dx = numpy.fft.rfft(d, nx)
450 f = real_fft_freq(dx.shape[0])
451 q = 1.0/(2*cutoff_freq**2)
452 ff = numpy.exp(-numpy.square(f) * q)
453 return numpy.fft.irfft(dx*ff, nx)[:n]
454
455 lowpass_sym_gaussian = lowpass_sym_Gaussian
456
457
458
460 """A time-symmetric filter with a magnitude response that
461 is the same as the Butterworth filter.
462 Cutoff freq is measured in cycles per point."""
463 assert len(d.shape) == 1
464 if cutoff_freq<=0 or cutoff_freq>10.0:
465 die.warn('Silly value of cutoff_freq: %g cycles per point.' % cutoff_freq)
466 n = d.shape[0]
467 nx = _filter_win_size(n, cutoff_freq)
468 dx = numpy.fft.rfft(d, nx)
469 f = real_fft_freq(dx.shape[0])
470 ff1 = (f/cutoff_freq)**(2*order)
471 ff = numpy.sqrt( ff1 / (1 + ff1) )
472 return numpy.fft.irfft(dx*ff, nx)[:n]
473
474
476 q = numpy.cos(numpy.arange(10000)*0.003*2*math.pi)
477 sq = numpy.sum(q*q)
478 q1 = lowpass_sym_butterworth(q, 0.1)
479 sq1 = numpy.sum(q1*q1)
480
481 assert 0.99*sq < sq1 < 1.00001*sq
482 q2 = lowpass_sym_butterworth(q, 0.05)
483 sq2 = numpy.sum(q2*q2)
484
485 assert 0.99*sq < sq2 < 1.00001*sq
486 q4 = lowpass_sym_butterworth(q, 0.003)
487 sq4 = numpy.sum(q4*q4)
488
489 assert 0.4*sq < sq4 < 0.6*sq
490 q5 = lowpass_sym_butterworth(q, 0.001)
491 sq5 = numpy.sum(q5*q5)
492
493 assert 0 < sq5 < 0.2*sq
494
495
497 q = numpy.cos(numpy.arange(10000)*0.003*2*math.pi)
498 sq = numpy.sum(q*q)
499 q1 = hipass_sym_butterworth(q, 0.0003)
500 sq1 = numpy.sum(q1*q1)
501
502 assert 0.99*sq < sq1 < 1.00001*sq
503 q4 = hipass_sym_butterworth(q, 0.003)
504 sq4 = numpy.sum(q4*q4)
505
506 assert 0.4*sq < sq4 < 0.6*sq
507 q5 = hipass_sym_butterworth(q, 0.03)
508 sq5 = numpy.sum(q5*q5)
509
510 assert 0 < sq5 < 0.2*sq
511
512
514 """Create pink noise where the power spectral density
515 goes as 1/f (except right at f=0).
516 """
517 m = near_win_size_real(near=2*n+2, tol=n//3+1)
518
519 rtmp = numpy.random.standard_normal(m)
520 itmp = numpy.random.standard_normal(m)
521 tmp = rtmp + 1j*itmp
522 f = real_fft_freq(m)
523 assert f[0] == 0.0
524 tmp[0] = 0.0
525 f[0] = 1.0
526 numpy.multiply(tmp, 1.0/numpy.sqrt( f ), tmp)
527 return numpy.fft.inverse_real_fft(tmp)[:n]
528
529 if __name__ == '__main__':
530 test_win_size()
531 test_print_win()
532 test_hpf()
533 test_lpb()
534 test_hpb()
535
536