forked from ShihPingLai/Group-3
-
Notifications
You must be signed in to change notification settings - Fork 0
/
produce_sig.py
91 lines (71 loc) · 1.63 KB
/
produce_sig.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
import numpy as np
import scipy.integrate as integrate
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
#define universal variables
c0 = 15.6
c1 = 1.0
c2 = 28.0
m0 = -1.143
m1 = -0.714
#open the data from our sound
f1= open("sound1.txt","r")
f2= open("sound2.txt","r")
f3= open("time.txt","r")
#just a little extra, quite unimportant
def f(x):
f = m1*x+(m0-m1)/2.0*(abs(x+1.0)-abs(x-1.0))
return f
#the actual function calculating
def dH_dt(H, t=0):
return np.array([c0*(H[1]-H[0]-f(H[0])),
c1*(H[0]-H[1]+H[2]),
-c2*H[1]])
#computational time steps
t=[]
sound1=[]
sound2=[]
for i in range (720832):
t.append(float(f3.readline()))
sound1.append(float(f1.readline()))
sound2.append(float(f2.readline()))
#x, y, and z initial conditions
H0 = [0.7, 0.0, 0.0]
H, infodict = integrate.odeint(dH_dt, H0, t, full_output=True)
print(infodict['message'])
#plt.subplot(221)
fig1=plt.figure(1)
ax = fig1.add_subplot(111, projection='3d')
ax.plot(H[:,0]*5000, H[:,1]*5000, H[:,2]*5000)
plt.show()
#plt.subplot(222)
"""
fig2=plt.figure(2)
x = np.sin(t)
plt.plot(t, x)
plt.show()
"""
#plt.subplot(223)
fig3 = plt.figure(3)
ax = fig3.add_subplot(111, projection='3d')
ax.plot(H[:,0]*5000+sound1, H[:,1]*5000+sound2, H[:,2]*5000)
plt.show()
out1=H[:,0]*5000+sound1
out2=H[:,1]*5000+sound2
#print (out)
fig4=plt.figure(4)
plt.plot(t,out1)
plt.show()
fig5=plt.figure(5)
plt.plot(t,sound2)
plt.show()
f4= open("in1.txt","w")
f5= open("in2.txt","w")
for i in range(720832):
f4.write(str(out1[i])+"\n")
f5.write(str(out2[i])+"\n")
f4.close()
f5.close()
f1.close()
f2.close()
f3.close()