-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathfft.py
executable file
·173 lines (166 loc) · 4.88 KB
/
fft.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
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
#!/usr/bin/env python
#coding=utf8
#========================================
# FFT and Velocity visualisation script
# Created by Max Dukov
# maxdukov@gmail.com
#========================================
# v. 1.0.090417 =)
print "FFT for LIS3DH vis script v1.0.140817"
import numpy as np
from numpy import array, arange, abs as np_abs
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from datetime import datetime
import sqlite3
import argparse
import socket
import fcntl
import struct
def get_ip_address(ifname):
s = socket.socket(socket.AF_INET, socket.SOCK_DGRAM)
return socket.inet_ntoa(fcntl.ioctl(
s.fileno(),
0x8915, # SIOCGIFADDR
struct.pack('256s', ifname[:15])
)[20:24])
freq = 400/2
ampl = 140000
#==============================
# setup arguments parsing here
#==============================
#parser = argparse.ArgumentParser()
#parser.add_argument("--sensor", type=int, help="sensor id")
#parser.add_argument("--norm", type=int, help="select normalization mode. 1 for normalization, 0 for standart(default)")
#args = parser.parse_args()
#if args.sensor is None:
# print "Default sensor id 0x1d will be used. Use ./fft.py --sensor 2 for select second sensor. Use ./fft.py -h for help."
# args.sensor = 1
#if (args.norm is None or args.norm == 0):
# print "Default mode, no normalization"
# norm = 0
#if args.norm == 1:
# print "Normalization mode"
# norm = 1
#=========================
# init sqlite connection
#=========================
con = sqlite3.connect('/var/www/html/gyro.db')
con.row_factory = lambda cursor, row: row[0]
cur = con.cursor()
cur.execute('SELECT x FROM log_lis ORDER BY dt ASC;')
x_ = cur.fetchall()
cur.execute('SELECT y FROM log_lis ORDER BY dt ASC;')
y_ = cur.fetchall()
cur.execute('SELECT z FROM log_lis ORDER BY dt ASC;')
z_ = cur.fetchall()
con.close()
print "==================================="
print "=> Processing data for LIS3DH sensor"
print "==================================="
######## process X
# freq graph
Yx = np_abs(np.fft.rfft(x_))/ampl
Yx_h = np_abs(np.fft.rfft(x_*np.hanning(len(x_))))/ampl
Yx[0] = 0
end = len(Yx)
X = np.asarray(np.linspace(0, freq, end, endpoint=True))
id = np.where( Yx == max(Yx))
print "max Yx=", max(Yx)
print "id=", id[0]
freq_x = X[id[0]] #
print "freq=", freq_x
print "=============="
### set plt size
fig_size = plt.rcParams["figure.figsize"]
fig_size[0] = 16
fig_size[1] = 22
plt.rcParams["figure.figsize"] = fig_size
fig = plt.figure()
fig.suptitle('FFT and velocity graphs for LIS3DH sensor ', fontsize=14, fontweight='bold')
### plot freq for X
plt.subplot(3, 2, 1)
plt.grid(True)
plt.plot(X,Yx,'b.-',label='Row FFT')
plt.plot(X,Yx_h,'c.-',label='FFT with Haning')
plt.xlabel('Freq (Hz)')
plt.legend(loc=2)
# Velocity graph
Vx = [0]*freq
for i in range(0,freq):
Vx[i] = x_[i]*(39.3701/(6.28*freq_x))
end = len(Vx)
X = np.asarray(np.linspace(0, 0.5, end, endpoint=True))
### plot Velocity for X
plt.subplot(3, 2, 2)
plt.grid(True)
plt.plot(X,Vx)
plt.xlabel('time (s)')
######## process Y
# freq graph
Yy = np_abs(np.fft.rfft(y_))/ampl
Yy_h = np_abs(np.fft.rfft(y_*np.hanning(len(y_))))/ampl
Yy[0] = 0
end = len(Yy)
Xy = np.asarray(np.linspace(0, freq, end, endpoint=True))
id = np.where( Yy == max(Yy))
print "max Yy=", max(Yy)
print "id=", id[0]
freq_y = Xy[id[0]] #
print "freq=", freq_y
print "=============="
### plot freq for Y
plt.subplot(3, 2, 3)
plt.grid(True)
plt.plot(Xy,Yy,'b.-',label='Row FFT')
plt.plot(Xy,Yy_h,'c.-',label='FFT with Hanning')
plt.legend(loc=2)
plt.xlabel('Freq (Hz)')
# Velocity graph
Vy = [0]*freq
for i in range(0,freq):
Vy[i] = y_[i]*(39.3701/(6.28*freq_y))
end = len(Vy)
Xv = np.asarray(np.linspace(0, 0.5, end, endpoint=True))
### plot Velocity for Y
plt.subplot(3, 2, 4)
plt.grid(True)
plt.plot(Xv,Vy)
plt.xlabel('time (s)')
######## process Z
# freq graph
Yz = np_abs(np.fft.rfft(z_))/ampl
Yz_h = np_abs(np.fft.rfft(z_*np.hanning(len(z_))))/ampl
Yz[0] = 0
end = len(Yz)
Xz = np.asarray(np.linspace(0, freq, end, endpoint=True))
id = np.where( Yz == max(Yz))
print "max Yz=", max(Yz)
print "id=", id[0]
freq_z = Xz[id[0]] #
print "freq=", freq_z
### plot freq for Z
plt.subplot(3, 2, 5)
plt.grid(True)
plt.plot(Xz,Yz,'b.-',label='Row FFT')
plt.plot(Xz,Yz_h,'c.-',label='FFT with Hanning')
plt.legend(loc=2)
plt.xlabel('Freq (Hz)')
# Velocity graph
Vz = [0]*freq
for i in range(0,freq):
Vz[i] = z_[i]*(39.3701/(6.28*freq_z))
end = len(Vz)
Xv = np.asarray(np.linspace(0, 0.5, end, endpoint=True))
### plot Velocity for Z
plt.subplot(3, 2, 6)
plt.grid(True)
plt.plot(Xv,Vz)
plt.xlabel('time (s)')
timestr = datetime.strftime(datetime.now(), '%Y-%m-%d_%H:%M:%S')
#if norm == 0:
plt.savefig('/var/www/html/fft_all_LIS_'+timestr+'_.png')
#if norm == 1:
# plt.savefig('/var/www/html/fft_all_norm_'+timestr+'_'+(str(args.sensor))+'.png')
print 'Please have a look on results here: http://'+(str(get_ip_address('eth0')))+'/fft_all_LIS_'+timestr+'.png'