import sys
from math import pi, sin, cos, tan, exp
import random
import numpy as np
from matplotlib import pyplot as plt
#x0, x1, nData = 0, 1.0, 128
x0, x1, nData = 0, 1.0, 1024
WG = 0.02
IgnoreZero = 0
ConvFunc = 'Gaussian'
mmx = 38
lmax = 50
[ドキュメント]
def MEMSub(nmax, mmx, dt, pg, pf, pspec, pm):
f0 = 1.0 / dt / (nmax - 1.0)
for i in range(1, nmax+1):
pf[i] = f0 * i
sumr = 1.0
sumi = 0.0
cr = 0.0
ci = 2.0 * pi * pf[i] * dt
for j in range(1, mmx+1):
sumr += pg[j] * cos(j * ci)
sumi += pg[j] * sin(j * ci)
pspec[i] = 2.0 * pm * dt / (sumr * sumr + sumi * sumi)
[ドキュメント]
def MEMLevins(nmax, pg, pm, pb1, pb2, pa, m):
sn, sd = 0.0, 0.0
for i in range(1, nmax+1):
sn += pb1[i] * pb2[i]
sd += pb1[i] * pb1[i] + pb2[i] * pb2[i]
pg[m] = -2.0 * sn / sd
pm *= (1.0 - pg[m] * pg[m])
if m != 1:
for k in range(1, m):
pg[k] = pa[k] + pg[m] * pa[m-k]
for i in range(1, m):
pa[i] = pg[i]
for i in range(1, nmax - m + 1):
pb1[i] += pa[m] * pb2[i]
pb2[i] = pb2[i+1] + pa[m] * pb1[i+1]
[ドキュメント]
def MEMBurg(nmax, mmx, lmax, px, pg, pc, pfpe):
nmat = max([mmx+2, nmax+2])
b1 = np.empty(nmat)
b2 = np.empty(nmat)
a = np.empty(nmat)
sum = 0.0;
for i in range(1, nmax + 1):
sum += px[i] * px[i]
pc[1] = sum / nmax
pm = pc[1]
pfpe[1] = (nmax + 1.0) / (nmax - 1.0) * pm
lag = mmx + 1
lg1 = lag + 1
b1[1] = px[1]
for i in range(2, nmax + 1):
b1[i] = px[i]
b2[i-1] = px[i]
for m in range(1, mmx + 1):
MEMLevins(nmax, pg, pm, b1, b2, a, m)
sum = 0.0;
for i in range(1, m + 1):
sum -= pc[m-i] * pg[i]
pc[m+1] = sum
if m != nmax - 1:
pfpe[m+1] = (nmax + m + 1) / (nmax - m - 1) * pm
if lag < lmax:
for l in range(lg1, lmax + 1):
sum = 0.0
for i in range(1, mmx + 1):
sum -= pc[l] * pg[i]
pc[l] = sum
return pm
[ドキュメント]
def MEM(dt, px, mmx, lmax):
nmax = len(px)
nmat = max([mmx+2, nmax+2])
f = np.zeros(nmat)
fpe = np.zeros(nmat)
spec = np.zeros(nmat)
c = np.zeros(nmat)
g = np.zeros(nmat)
x2 = [0.0] + px
pm = MEMBurg(nmax, mmx, lmax, x2, g, c, fpe)
MEMSub(nmax, mmx, dt, g, f, spec, pm)
return f, fpe, spec, c, mmx, lmax, pm
[ドキュメント]
def MEMCall(x, y, mmx, lmax):
lag = mmx + 1
if lag+1 > lmax:
lmax = lag + 1
nmax = len(x)
sum = 0.0;
for i in range(nmax):
sum += x[i]
av = sum / nmax
for i in range(nmax):
x[i] -= av
dt = x[1] - x[0]
pf, pfpe, pspec, pc, mmx , lmax, pm = MEM(dt, y, mmx , lmax)
lag = mmx + 1
lg1 = lag + 1
lx1 = lmax + 1
print("\nf\tspc\tfpe\tauto c.")
for i in range(lag):
print("{}\t{}\t{}\t{}".format(pf[i] , pspec[i] , pfpe[i] , pc[i]))
if lag != lmax:
for i in range(lg1, lmax+1):
print("{}\t{}\t{}".format(pf[i], pspec[i], pc[i]))
for i in range(lx1, nmax):
print("{}\t{}".format(pf[i], pspec[i]))
return pf, pspec, pfpe, pc
[ドキュメント]
def Func(x):
x += 0.03 * random.random()
f1, p1, A1 = 1.5, pi/4.0, 1.0
f2, p2, A2 = 3.0, pi/3.0, 0.3
f3, p3, A3 = 10.0, pi/6.0, 0.5
return A1 * sin(2.0*pi * f1 * x + p1) \
+ A2 * sin(2.0*pi * f2 * x + p2) \
+ A3 * sin(2.0*pi * f3 * x + p3)
[ドキュメント]
def main():
xstep = (x1 - x0) / nData
print("x = ({}, {}, {})".format(x0, x1, xstep))
print("nData=", nData)
print("W(Gauss)=", WG)
x = []
y = []
for i in range(nData):
xx = x0 + i * xstep
x.append(xx)
y.append(Func(xx))
pf, pspec, pfpe, pc = MEMCall(x, y, mmx, lmax)
#=============================
# Plot graphs
#=============================
fig = plt.figure()
ax1 = fig.add_subplot(2, 1, 1)
ax2 = fig.add_subplot(2, 1, 2)
ax1.plot(x, y, label = 'input')
ax1.set_xlabel("t")
ax1.set_ylabel("y")
ax1.legend()
ax2.plot(pf, pspec, label = 'mem')
ax2.set_xlabel("f")
ax2.set_ylabel("mem")
ax2.legend()
plt.tight_layout()
plt.pause(0.1)
print("Press ENTER to exit>>", end = '')
input()
if __name__ == '__main__':
main()