00001 """Compute maximum entropy coefficients for data. Based loosely on the the
00002 concepts in "Numerical Recipes in C", 2nd ed., by Press, Flannery, Teukolsky and Vetterling (q.v.)
00003 but I don't think it is any copyright infringement"""
00004 _rcsid="$Id: memcof.py,v 1.3 2003/05/30 13:31:55 mendenhall Release-20030716 $"
00005
00006 import Numeric
00007 import math
00008
00009 def memcof(data, poles):
00010 n=len(data)
00011 xms=Numeric.dot(data,data)/float(n)
00012 wk1=data[:-1]
00013 wk2=data[1:]
00014 d=Numeric.zeros(poles,Numeric.Float)
00015
00016 for k0 in range(poles):
00017 num=Numeric.dot(wk1,wk2)
00018 denom=Numeric.dot(wk1, wk1)+Numeric.dot(wk2, wk2)
00019 d[k0]=2.0*num/denom
00020 xms*=(1.0-d[k0]**2)
00021 if k0!=0:
00022 d[:k0]=wkm-d[k0]*wkm[-1::-1]
00023 if k0!=poles-1:
00024 wkm=d[:k0+1]
00025 wk1, wk2 = wk1[:-1]-wkm[k0]*wk2[:-1], wk2[1:]-wkm[k0]*wk1[1:]
00026
00027 return xms, d
00028
00029
00030 def evlmem(fdt, d, xms):
00031 n=len(d)
00032 theta=2*math.pi*fdt*(Numeric.array(range(1,n+1), Numeric.Float))
00033 zr=1.0-Numeric.dot(Numeric.cos(theta), d)
00034 zi=Numeric.dot(Numeric.sin(theta), d)
00035 return xms/(zr*zr+zi*zi)
00036
00037
00038 if __name__=="__main__":
00039
00040 datalen=10000
00041 poles=20
00042 zfreq=0.21345
00043 pspoints=200
00044 damping=-20.0/datalen
00045 noise=1.0
00046
00047 xvals=Numeric.array(range(datalen),Numeric.Float)
00048 data=Numeric.sin(xvals*(2.0*math.pi*zfreq)) * Numeric.exp(xvals*damping)
00049 import random
00050 r=random.Random(1234)
00051 data+= Numeric.array([r.random()-0.5 for i in range(datalen)])*noise
00052
00053 d2=Numeric.dot(data,data)/datalen
00054
00055 xms, d = memcof(data, poles)
00056
00057 freqlist = [0.5*i/pspoints for i in range(pspoints)]
00058 pspect = [evlmem(f, d, xms) for f in freqlist]
00059
00060 pssum=2.0*Numeric.sum(pspect)*(0.5/pspoints)
00061
00062 print "input power = ", d2, "output power = ", pssum, "ratio =", pssum/d2
00063 import graphite
00064 g=graphite.Graph()
00065 d=graphite.Dataset()
00066 d.x=freqlist
00067 d.y=pspect
00068 g.datasets=d
00069 graphite.genOutput(g,'QD',canvasname="Power Spectrum")
00070
