[frames] | no frames]

# Source Code for Module lib.voice_misc

```  1  #!/usr/bin/env python
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
12 -class cont_kernel_c(ortho_poly.ortho):
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          """
19 -        def __init__(self, n=None, x=None):
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
25 -        def wt_(self):
26                  return numpy.square(1.0-numpy.square(self.x))
27
28 -        def compute(self, nnodes):
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
41 -def window(n, k, norm=2):
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          # opn.x ranges over (-1, 1)
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
63 -def test_cont_kernel():
64          import pylab
65          c = cont_kernel_c(20)
66          pylab.plot(c.P(0))
67          pylab.plot(c.P(1))
68          pylab.plot(c.P(2))
69          pylab.plot(c.P(3))
70          pylab.show()
71
72
73 -def cont_kernel(n, k, norm=2, kernel=window):
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)
97          if norm > 0:
98                  numpy.divide(w, (numpy.absolute(w)**norm).sum()**(1.0/norm), w)
99          return w
100
101
102 -def xexp(a):
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
111 -def start_end(t, dt, window, len_d):
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
134 -def _filter_win_size(n, cutoff_freq):
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 = {}
142 -def _fftwork(n):
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
161 -def _fftwork_real(n):
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                          # This branch is just to keep near_win_size
204                          # from trying too hard.   Otherwise, it's quite
205                          # possible to spend more time finding the optimal
206                          # size than one spends taking the actual FFT.
207                          since += 1
208                          if since > Since_lim:
209                                  break
210          return besti
211
212
213 -def frommid_generator(s, e):
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
223 -def near_win_size_real(near=None, tol=None, min=None, max=None):
224          return near_win_size(near, tol, min, max, real=True)
225
226
227 -def test_win_size():
228          assert near_win_size(min=4, max=7) == 4
229          assert near_win_size(min=5, max=7) != 7
230          assert near_win_size(near=15, tol=2) != 17
231          assert near_win_size(near=14, tol=1) != 13
232          assert near_win_size(near=97, tol=1) != 97
233          assert near_win_size(near=197, tol=1) != 197
234          assert near_win_size(near=23, tol=1) != 23
235
236
237 -def test_print_win():
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  # def hipass_first_order(d, cutoff_freq):
251          # """A first-order high-pass filter.
252          # Cutoff freq is measured in cycles per point."""
253  #
254          # assert cutoff_freq >= 0.0
255          # assert len(d.shape) == 1
256          # n = d.shape[0]
257          # nx = _filter_win_size(n, cutoff_freq)
258          # dx = Num.FFT.rfft(d, nx)
259          # f = real_fft_freq(dx.shape[0])
260          # fr = 1j* (f / cutoff_freq)
261          # ff = fr / (1 + fr)
262          # ff[0] = ff[0].real
263          # ff[-1] = ff[-1].real
264          # return Num.FFT.irfft(dx*ff, nx)[:n]
265  #
266  #
267  # def lopass_first_order(d, cutoff_freq):
268          # """A first-order low-pass filter.
269          # Cutoff freq is measured in cycles per point."""
270  #
271          # assert cutoff_freq >= 0.0
272          # assert len(d.shape) == 1
273          # n = d.shape[0]
274          # nx = _filter_win_size(n, cutoff_freq)
275          # dx = Num.FFT.rfft(d, nx)
276          # f = real_fft_freq(dx.shape[0])
277          # fr = 1j* (f / cutoff_freq)
278          # ff = 1 / (1 + fr)
279          # ff[0] = ff[0].real
280          # ff[-1] = ff[-1].real
281          # return Num.FFT.irfft(dx*ff, nx)[:n]
282
283  hipass_first_order = GT.hipass_first_order
284  lopass_first_order = GT.lopass_first_order
285
286
287 -def bandpass_first_order(d, cutoff_hp, cutoff_lp):
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
309 -def test_hpf():
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          # print sq5/sq
321          assert  0.003 < sq5 < 0.03*sq
322
323
324
325
326 -def fft_freq(n, d=1.0):
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
340 -def fft_indices(f0, f1, n, d=1.0):
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
347 -def real_fft_freq(n, d=1.0):
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
362 -def real_fft_indices(f0, f1, n, d=1.0):
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
374 -def dataprep_flat_real(data, dt, edgewidth, pad=None):
375          return dataprep_flat_generic(data, dt, edgewidth, pad, True)
376
377 -def dataprep_flat(data, dt, edgewidth, pad=None):
378          return dataprep_flat_generic(data, dt, edgewidth, pad, False)
379
380
381 -def dataprep_flat_generic(data, dt, edgewidth, pad, isreal):
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]
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
402 -def dataprep2_flat_generic(data, dt, edgewidth, pad, isreal):
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]
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
423 -def lowpass_sym_butterworth(d, cutoff_freq, order=4):
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
440 -def lowpass_sym_Gaussian(d, cutoff_freq):
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
459 -def hipass_sym_butterworth(d, cutoff_freq, order=4):
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
475 -def test_lpb():
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          # print sq1/sq
481          assert 0.99*sq < sq1 < 1.00001*sq
482          q2 = lowpass_sym_butterworth(q, 0.05)
483          sq2 = numpy.sum(q2*q2)
484          # print sq2/sq
485          assert 0.99*sq < sq2 < 1.00001*sq
486          q4 = lowpass_sym_butterworth(q, 0.003)
487          sq4 = numpy.sum(q4*q4)
488          # print sq4/sq
489          assert 0.4*sq < sq4 < 0.6*sq
490          q5 = lowpass_sym_butterworth(q, 0.001)
491          sq5 = numpy.sum(q5*q5)
492          # print sq5/sq
493          assert  0 < sq5 < 0.2*sq
494
495
496 -def test_hpb():
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          # print sq1/sq
502          assert 0.99*sq < sq1 < 1.00001*sq
503          q4 = hipass_sym_butterworth(q, 0.003)
504          sq4 = numpy.sum(q4*q4)
505          # print sq4/sq
506          assert 0.4*sq < sq4 < 0.6*sq
507          q5 = hipass_sym_butterworth(q, 0.03)
508          sq5 = numpy.sum(q5*q5)
509          # print sq5/sq
510          assert  0 < sq5 < 0.2*sq
511
512
513 -def pink_noise(n):
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          # test_cont_kernel()
536
<!--
expandto(location.href);
// -->

```

 Generated by Epydoc 3.0.1 on Thu Jun 16 20:03:15 2011 http://epydoc.sourceforge.net