Commits

Jason S committed 447e173

updated to fix minor bugs, added sigmadeltatest

  • Participants
  • Parent commits 1b844ce

Comments (0)

Files changed (2)

sigmadeltatest.py

+#   Copyright 2013 Jason Sachs
+#
+#   Licensed under the Apache License, Version 2.0 (the "License");
+#   you may not use this file except in compliance with the License.
+#   You may obtain a copy of the License at
+#
+#       http://www.apache.org/licenses/LICENSE-2.0
+#
+#   Unless required by applicable law or agreed to in writing, software
+#   distributed under the License is distributed on an "AS IS" BASIS,
+#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#   See the License for the specific language governing permissions and
+#   limitations under the License.
+
+from twopyf import *
+import numpy
+import matplotlib.pyplot as plt
+
+def sigmadeltamod1(xit, i1_0=0):
+  x2 = (2*(x-0.5) for x in xit)
+  feedback = Deferred()
+  integ = integrator(subtract(x2,feedback), i1_0, immediate=True)
+  (integ, integout) = itertools.tee(integ)
+  comp = comparator(integ, 0)
+  (comp, compout) = itertools.tee(comp)
+  feedback.set(unitdelay(switch(comp,-1,1)))
+  return (compout, integout)
+
+def sigmadeltamod2(xit, k, i1_0=0, i2_0=0):
+  x2 = (2*(x-0.5) for x in xit)
+  feedback = Deferred()
+  (feedback0, feedback1) = itertools.tee(buffer(feedback))
+  i1 = integrator(scale(subtract(x2,feedback0), k), i1_0, immediate=True, limits=[-1,1])
+  (i1, i1out) = itertools.tee(i1)
+  i2 = integrator(subtract(i1,feedback1), i2_0, immediate=True)
+  (i2, i2out) = itertools.tee(i2)
+  (y,yfb) = itertools.tee(comparator(i2,0))
+  feedback.set(unitdelay(switch(yfb,-1,1)))
+  return (y, i1out, i2out)
+
+def atruncater(n):
+  def f(x):
+    return numpy.array(list(truncate(x,n)))
+  return f
+
+def filter1pk(x,alpha,k):
+  y = x
+  for i in range(k):
+    y = filter1p(y,alpha)
+  return y
+
+dt = 100e-6
+t = numpy.arange(0,1,dt)
+Tchop = 1.0/11
+n = len(t)
+x = ((t*(numpy.mod(t/Tchop,2)<1))-0.5)*0.9+0.5
+
+(d,di1) = sigmadeltamod1(x)
+k2 = 0.1
+(d2,di2a,di2b) = sigmadeltamod2(x,k2)
+snip = atruncater(n) 
+tau = 0.005
+alpha = dt/tau
+#xfilt1 = snip(filter1pk(x),alpha,2))
+xfilt1 = snip(filter1pk(x,alpha,2))
+di1 = snip(di1)
+y1 = snip(filter1pk(d,alpha,2))
+
+xfilt2 = snip(filter1pk(x,alpha,2))
+di2a = snip(di2a)
+di2b = snip(di2b)
+y2 = snip(filter1pk(d2,alpha,2))
+
+for (i, data) in enumerate([[xfilt1,y1,di1],[xfilt2,y2,di2a,di2b]]):
+  xfilt = data[0]
+  y = data[1]
+  n = len(data)
+  fig = plt.figure(i)
+  ax=fig.add_subplot(n,1,1)
+  ax.plot(t,x,t,xfilt,t,y)
+  for j in range(n-2):
+    ax=fig.add_subplot(n,1,j+2)
+    ax.plot(t,data[j+2])
+  ax=fig.add_subplot(n,1,n)
+  ax.plot(t,xfilt-y)
+  break
+plt.show()
 import itertools
 import random
 
-def buffer(x):
-  for item in x:
-    yield item
+def trace(g, fmt):
+  for x in g:
+    print fmt % x
+    yield x
 
-def scale(x, k):
-  for item in x:
-    yield k*item
+def buffer(xit):
+  for x in xit:
+    yield x
 
-def unitdelay(x, x0):
+def bufferF(xit, f):
+  for x in xit:
+    yield f(x)
+
+def scale(xit, k):
+  for x in xit:
+    yield k*x
+
+def unitdelay(xit, x0=0):
   yield x0
-  for item in x:
-    yield item
+  for x in xit:
+    yield x
+
+def merge(g1, g2):
+  try:
+    i1 = iter(g1)
+    i2 = iter(g2)
+    while True:
+      yield (i1.next(), i2.next())
+  except StopIteration:
+    pass
 
 def add(g1, g2):
-  for (item1, item2) in itertools.izip(g1,g2):
-    yield item1+item2
+  for (x1, x2) in itertools.izip(g1,g2):
+    yield x1+x2
 
 def subtract(g1, g2):
-  for (item1, item2) in itertools.izip(g1,g2):
-    yield item1-item2
+  for (x1, x2) in itertools.izip(g1,g2):
+    yield x1-x2
     
-def gaussnoise(mu, sigma, seed=None):
+def infiniteStream(x0=0):
+  while True:
+    yield x0
+
+def truncate(xit, n):
+  for x in xit:
+    if n <= 0:
+      break
+    n -= 1
+    yield x
+    
+def gaussnoise(mu=0, sigma=1, n=None, seed=None):
   r = random.Random()
   r.seed(seed)
-  while True:
-    yield r.gauss(mu, sigma)    
+  it = infiniteStream() if n is None else xrange(n)
+  for x in it:
+    yield r.gauss(mu, sigma)       
     
-def integrator(x, x0=0):
+def integrator(xit, x0, immediate=False, limits=None):
   S = x0
-  yield S
-  for item in x:
-    S += item
+  if limits is None:
+    ulim = float('inf')
+    llim = -ulim
+  else:
+    (llim, ulim) = limits
+  if not immediate:
+    yield S
+  for x in xit:
+    S += x
+    S = min(ulim,max(llim,S))
     yield S
   
-def differentiator(x, x0=0):
+def differentiator(xit, x0):
   xprev = x0
-  for item in x:
-    yield item - xprev
-    xprev = item
+  for x in xit:
+    yield x - xprev
+    xprev = x
+
+def comparator(xit, threshold=0):
+  for x in xit:
+    yield 1 if x >= threshold else 0
+    
+def switch(xit, *args):
+  for x in xit:
+    yield args[x]
+    
+def filter1p(xit, alpha, y0=0):
+  y = y0
+  yield y
+  for x in xit:
+    y += alpha*(x-y)
+    yield y
+
+class Deferred(object):
+  def __init__(self, f=None, *args):
+    self.f = f if f is not None else iter
+    self.args = args
+  def set(self, value):
+    self.value = value
+  def __iter__(self):
+    return self.f(self.value, *self.args)
   
-def comparator(x, threshold):
-  for item in x:
-    yield 1 if item >= threshold else 0
+#####
+
+if __name__ == '__main__':
+  x = xrange(10)
+  f = Deferred()
+  u = unitdelay(f)
+  [u,u2] = itertools.tee(u)
+  f.set(add(u2,x))
+  for item in u:
+    print "item=%s" % item