Source Code for Module gmisclib.sharp_energy

 1  import Num 
 2  import math 
 3  import cmath 
 4  import dpss 
 5  import gpkmisc 
 6   
 7   
 8   
 9   
10 -def complex_median(x): 
11          if x.typecode() == Num.Float: 
12                  return gpkmisc.N_median(x) 
13   
14          r1 = gpkmisc.N_median(x.real) 
15          i1 = gpkmisc.N_median(x.imag) 
16          phasor = cmath.exp(1j*math.pi/4) 
17          t = x*phasor 
18          r2 = gpkmisc.N_median(t.real) 
19          i2 = gpkmisc.N_median(t.imag) 
20          return 0.5*(complex(r1,i1) + complex(r2,i2)/phasor) 
21   
22   
23 -def one_median(binw, center, v, nbins=None): 
24          """Local power in a signal. 
25          dt and binw are measured in bins. 
26          """ 
27          if nbins is None: 
28                  nbins = 3 
29          assert type(binw) == type(1) 
30          assert binw>=1 
31          hbw = binw/2.0 
32   
33          if center == "?data length":    # How much data do you need? 
34                  return int(math.ceil(hbw*(4*nbins+1))) + binw 
35   
36          w = 0.5/binw 
37          # print "binw=", binw, "1/w=", 1/w, "nbins=", nbins, "center=", center 
38          lmbda, window = dpss.dpsscache(w, 1, binw) 
39          assert lmbda[0]>0.7  and  lmbda[0]<=1.0 
40          v = Num.asarray(v) 
41          t = Num.zeros((4*nbins+1,), v.typecode()) 
42          m = 0 
43          for k in range(-2*nbins, 2*nbins+1): 
44                  ctr = center + hbw * k 
45                  s = int(round(ctr - hbw)) 
46                  if s < 0: 
47                          continue 
48                  e = s + binw 
49                  if e >= len(v): 
50                          continue 
51                  t[m] = Num.sum(v[s:e] * window[0]) 
52                  # print "bin", k, "s=", s, "e=", e, "t=", t[m] 
53                  m += 1 
54          return complex_median(t[:m]) / Num.sum(window[0]) 
55   
56   
57   
58 -def stepped_median(dt, binw, v, nbins=None): 
59          """Local power in a signal. 
60          dt and binw are measured in bins. 
61          """ 
62          binw = int(round(binw)) 
63          v = Num.asarray(v) 
64          nb = int(round(len(v)/dt)) 
65          o = Num.zeros((nb,), v.typecode()) 
66          for i in range(nb): 
67                  o[i] = one_median(binw, dt*i, v, nbins) 
68          return o 
69   
70   
71 -def pwr(dt, binw, v, nbins=3): 
72          v = Num.asarray(v, Num.Float) 
73          return stepped_median(dt, binw, Num.square(v), nbins) 
74   
75   
76  if __name__ == '__main__': 
77          print pwr(1, 3, [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]) 
78          # print pwr(1, 3, [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]) 
79          # print pwr(1, 3, [0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1]) 
80