[frames] | no frames]

# 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
<!--
expandto(location.href);
// -->

```

 Generated by Epydoc 3.0.1 on Thu Sep 22 04:25:13 2011 http://epydoc.sourceforge.net