This repository was archived by the owner on Feb 28, 2022. It is now read-only.
-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathAnalysis_3.py
More file actions
312 lines (257 loc) · 9.56 KB
/
Analysis_3.py
File metadata and controls
312 lines (257 loc) · 9.56 KB
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
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
#!/usr/bin/python3.4
##
# Analysis.py
#
# Author: Vincent Steffens, vsteffen@ucsc.edu
# Date: 16 November 2014
#
# Produces a mean power spectrum of raw SEAD plug current
# data, normalized by total current.
# Outputs to spectrum in a numpy array to stdout or a text
# file, each array element on a line.
##
#For numerical analysis
import numpy as np
import math
#For visualization
import matplotlib.pyplot as pyp
#For manipulating command-line arguments
import sys
#For handling files
import os
#For using regular expressions
import re
blockwidth = 200
Times = []
def init():
#Check for proper number of arguments, die if necessary
#max args = 2 + number of valid options
max_args = 10
if len(sys.argv) < 2:
pass
#print "Usage: Analysis.py [-cfhpvVw] source_file"
#print "Call with -h for help"
#sys.exit(-1)
if len(sys.argv) > max_args:
pass
#print "Usage: Analysis.py [-cfhpvVw] source_file"
#print "Call with -h for help"
#sys.exit(-2)
#handle calling with only Analysis.py -h
#if sys.argv[1] == '-h':
pass
#print_help()
#exit()
#Check for options. Break the option area of argv into a list of
#characters. Make a list of those that are valid, and a list of
#those that aren't. Let the user know if they endered an invalid
#option.
#When adding options, be sure to add a description to the -h section
#below
valid_options = 'cdfhpvVw'
alleged_options = list(set(''.join(sys.argv[1:-1])))
options = [x for x in alleged_options if re.search(x, valid_options)]
non_options = [x for x in alleged_options if x not in options and x != '-']
for i in non_options:
print("Ignoring invalid option: \'%s\'" % (i))
return options
#consider making times and currents an associative array
def import_and_trim():
Currents = []
amp_ids = [70, 66, 74, 62, 78, 58, 50, 14, 54]
#Try to open source file for reading
filename = sys.argv[len(sys.argv) - 1]
if os.path.isfile(filename):
with open(filename) as f:
#Check the first element of the first line.
#if it's "sensor_id", it's the old format
#if it's "1", it's the new format
line = f.readline().split(',')
#New format
if len(line) == 2:
#discard first one
for line in f:
line = line.split(',')
Currents.append(line[1])
elif line[0] == '1':
if line[1] == 'I':
Currents.append(line[3])
for line in f:
line = line.split(',')
if line[1] == 'I':
Currents.append(line[2])
else:
for line in f:
line = re.split(',|\t', line)
if int(line[0]) in amp_ids:
Currents.append(line[1])
else:
print("Analysis: cannot open file:", filename)
#Convert time since Unix epoch to intervals in microseconds
#Convert currents to milliamps
# Times = [ int(x) - int(Times[0]) for x in Times ]
Currents = [ 27*float(x)/1000 for x in Currents ]
return Currents
def produce_blocklist(Currents):
global blockwidth
blocklist = []
data_length = len(Currents) #or Times, just to de-specify
i = 0
while i < data_length:
if i + 200 > len(Currents):
break
blocklist.append(Currents[i:i+blockwidth])
i += blockwidth
return blocklist
def produce_mean_normalized_power_spectrum(blocklist):
#question: mean then normalize, or normalized then mean?
#doesn't matter because multiplication commutes
#units are in square amps
#here's why:
#original was in amps
#take ft (not fft)
#break up original into frequency components
#view from frequency perspective
#frequency components are in amps
#take mod square using absolute
#mod square(a + bi) = a^2 + b^2
#units are now in square amps per hz
#integrate it in piecese from a to a + blockwidth
#now you have units of amps sqared
#put those in an array
#oh yeah
#so now we have a power spectrum for an amp signal (odd sounding)
#in amps squared
#Produce power spectrum, sum
power_spectrum = np.square(np.absolute(np.fft.rfft(blocklist[0])))
sum_power_spectrum = power_spectrum
for i in range(1, len(blocklist)):
power_spectrum = np.square(np.absolute(np.fft.rfft(blocklist[i])))
sum_power_spectrum = np.add(sum_power_spectrum, power_spectrum)
#Take integral
total_sq_amps = 0.0
for i in range(0, len(sum_power_spectrum)):
total_sq_amps += sum_power_spectrum[i]
#Normalize by total amps measured
normalized_power_spectrum = sum_power_spectrum * (1/total_sq_amps)
#Finish taking mean
mean_normalized_spectrum = normalized_power_spectrum / len(blocklist)
return mean_normalized_spectrum
def display(spectrum):
template = np.ones(len(spectrum))
#Get the plot ready and label the axes
pyp.plot(spectrum)
max_range = int(math.ceil(np.amax(spectrum) / standard_deviation))
for i in range(0, max_range):
pyp.plot(template * (mean + i * standard_deviation))
# pyp.xlabel('Units?')
pyp.ylabel('Amps Squared')
pyp.title('Mean Normalized Power Spectrum')
if 'V' in Options:
pyp.show()
if 'v' in Options:
tokens = sys.argv[-1].split('.')
filename = tokens[0] + ".png"
input = ''
if os.path.isfile(filename):
input = input("Error: Plot file already exists! Overwrite? (y/n)\n")
while input != 'y' and input != 'n':
input = input("Please enter either \'y\' or \'n\'.\n")
if input == 'y':
pyp.savefig(filename)
else:
print("Plot not written.")
else:
pyp.savefig(filename)
def write_output():
tokens = sys.argv[-1].split('.')
filename = tokens[0] + ".txt"
#If a file with the same name already exists,
#check before overwriting and skip if necessary
if os.path.isfile(filename):
input = input("Error: Output file already exists! Overwrite? (y/n) : ")
while input != 'y' and input != 'n':
input = input("Please enter either \'y\' or \'n\'.\n")
if input == 'n':
print("Writing skipped.")
return
#Write
out = open(filename, 'w')
for element in Spectrum:
out.write(str(element) + ",")
out.close()
def print_help():
print("\nAnalysis.py.")
print("Used to produce a mean power spectrum of raw SEAD plug current data, normalized by total current")
print("\nCall using this syntax:")
print("$ python Analysis.py [-fhpvVw] [input file]")
print("")
print("Input: Raw, 4-column SEAD plug data")
print("Output: Either or both of:")
print("1. Numpy array printed to stdout.")
print("2. Spectrum written to text file.")
print("\nOptions may be arranged in any order, and in any number of groups")
print("")
print("-c:\t Classify. Match input data to a known group.")
print(" This only works for microwaves, lamps, and laptop computers")
print("-f:\t Fragmented data. This handles gaps in the data.")
print("-h:\t Help. Display this message.")
print("-p:\t Print. Print numpy array containing spectrum to terminal.")
print("-V:\t View. Display plot of spectrum, with the mean and multiples of the standard deviation.")
print("-v:\t Visualize. Print plot to file.")
print("-w:\t Write. Write spectrum to file, each array element on its own line")
print("\nExamples:")
print("1: Handle fragmented data, view plot, write spectrum to file")
print(" python Analysis.py -vfw 1_raw.csv")
print(" python Analysis.py -f -wv 1_raw.csv")
print("2: View plot of spectrum")
print(" python Analysis.py -V 1_raw.csv")
def classify(spectrum, mean, standard_deviation):
variation_coefficient = standard_deviation / mean
#count regions
count = 0
flag = 0
threshold = mean + standard_deviation
for i in range(0, len(spectrum)):
if flag == 0 and spectrum[i] > threshold:
flag = 1
count += 1
continue
elif flag == 1 and spectrum[i] < threshold:
flag = 0
output_string = sys.argv[-1]
ret = "Unknown Device"
if len(output_string) < 24:
output_string = output_string + "\t"
if variation_coefficient < 3.511:
#print "%s\t was recorded from a laptop computer." % (output_string)
ret = "computer"
elif np.amax(spectrum) < 0.01 or count > 1:
#print "%s\t was recorded from a microwave." % (output_string)
ret = "microwave"
else:
#print "%s\t was recorded from a lamp." % (output_string)
ret = "microwave"
return ret
def run(current_list):
Currents = []
Options = init()
Currents = current_list
Blocklist = produce_blocklist(Currents)
Spectrum = produce_mean_normalized_power_spectrum(Blocklist)
#mean and std are used by both display() and classify()
#only calculate once.
mean = np.mean(Spectrum)
standard_deviation = np.std(Spectrum)
return classify(Spectrum, mean, standard_deviation)
#This should be done first
if 'h' in Options:
print_help()
if 'c' in Options:
classify(Spectrum, mean, standard_deviation)
if 'p' in Options:
print(Spectrum)
if 'w' in Options:
write_output()
if 'v' in Options or 'V' in Options:
display(Spectrum)