-
Notifications
You must be signed in to change notification settings - Fork 1
/
sho3.py
110 lines (87 loc) · 2.66 KB
/
sho3.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
import random
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from scipy.integrate import odeint
from scipy.fft import fft, fftfreq
n = 10
k = 1000
b = 1
T = 1000
timestep = 0.01 # sec
window = 500 # sec
window_size = int(window/timestep)
def f(x, t, k, b):
ret = []
ret.append(0)
ret.append(0)
for i in range(n-2):
ret.append(x[2*i+3])
ret.append(k*x[2*i]-k*2*x[2*i+2]+k*x[2*i+4]+b *
(x[2*i]-x[2*i+2])**2+b*(x[2*i+4]-x[2*i+2])**2)
ret.append(0)
ret.append(0)
return ret
x0 = []
for i in range(n):
x0.append(i)
if i == 0 or i == n-1:
x0.append(0)
else:
x0.append(random.uniform(0, 1))
t = np.arange(0, T, step=timestep)
x = odeint(f, x0, t, args=(k, b,))
def do_fft(a=0, b=len(t)):
x_sample = x[a:b, :]
E_x = np.zeros(x_sample.shape[0])
freq = fftfreq(x_sample.shape[0], d=timestep)
for i in range(n):
y = np.abs(fft(x_sample[a:b, 2*i]-x0[2*i], n=x_sample.shape[0]))
E_x = E_x + (y**2)*(freq**2)*(2*np.pi)**2/2
return freq, E_x
fig = plt.figure(dpi=300)
ax = fig.gca()
#freq, E_x = do_fft(85000,90000)
freq_, E_x_ = do_fft(0, window_size)
# plt.plot(freq/freq[1],E_x,linestyle='-')
plt.plot(freq_/freq_[1], E_x_, linestyle='--')
#plt.show()
ax.set_xlim((np.min(freq_/freq_[1]), np.max(freq_/freq_[1])))
ax.set_ylim((0, np.max(E_x_)))
curve, = ax.plot([], [], lw=2, color='blue')
title = ax.text(0.5, 0.95, '', transform=ax.transAxes,
ha='center', va='center', fontsize=12)
def update(i):
freq, E_x = do_fft(i, i+window_size)
curve.set_data(freq/freq[1]*len(t)/window_size, E_x)
title.set_text('t='+str(i)+'~'+str(i+10)+'s')
return curve, title,
ani = animation.FuncAnimation(fig=fig, func=update, frames=len(
t)-window_size, interval=1, blit=True, repeat=True)
plt.show()
'''
fig = plt.figure(figsize=(4, 2), dpi=200)
ax = fig.gca()
ax.set_xlim((-0.5, n-0.5))
ax.axis('off')
title = ax.text(0.5, 0.65, '', transform=ax.transAxes,
ha='center', va='center', fontsize=12)
dot = []
for i in range(n):
d, = ax.plot([], [], color='black', marker='|' if i == 0 or i == n -
1 else 'o', markersize=15 if i == 0 or i == n-1 else 5, linestyle='')
dot.append(d)
def update(i):
title.set_text('t='+str(i))
for j in range(n):
dot[j].set_data(x[i][2*j], 0)
return tuple(dot)+(ax,)
def init():
title.set_text('t=0')
for j in range(n):
dot[j].set_data(x[0][2*j], 0)
return tuple(dot)+(ax,)
ani = animation.FuncAnimation(
fig=fig, func=update, frames=100000, init_func=init, interval=1, blit=True, repeat=True)
plt.show()
'''