Files
controller-wisetopapiserver/src/controllers/analysis/spikeDetect.py
T

110 lines
3.7 KiB
Python

import sys
import ast
import psycopg2
import numpy as np
import pandas as pd
from math import ceil, floor
from json import loads, dumps
from scipy.fft import fft, fftfreq, rfft, rfftfreq, irfft, ifft
if __name__ == '__main__':
params = {}
for input in sys.argv[1:]: #Now we're going to iterate over argv[1:] (argv[0] is the program name)
param = input.split("=") #Get what's left of the '='
params[param[0]] = param[1]
meta_id = params['meta_id']
channel = params['channel']
cut_off_freq = int(params['cut_off_freq'])
threshold = ast.literal_eval('[' + params['threshold'] + ']')
waveform = ast.literal_eval('[' + params['waveform'] + ']')
SAMPLE_RATE = int(params['sample_rate'])
DURATION = ceil(float(params['time_duration']))
DELTA_TIME = 1e6 / SAMPLE_RATE
N = SAMPLE_RATE * DURATION
ret_data_type = params['ret_data_type']
# print(meta_id, channel, cut_off_freq, threshold, waveForm, SAMPLE_RATE, DURATION, N)
conn = psycopg2.connect(database="postgres", user="biopro", password="BioProControlBox", host="127.0.0.1", port="5432")
cur = conn.cursor()
sql_str = 'SELECT data FROM "public"."' + str(channel) + '_recording_data_raws" WHERE parent = ' + meta_id + 'ORDER BY id ASC'
cur.execute(sql_str)
ret = None
try:
ret = cur.fetchall()
except BaseException as e:
print(e)
exit
finally:
time = np.array([], dtype=int)
data = np.array([])
for sub_data in ret:
for split_data in sub_data[0].split('"***"')[:-1]:
value = np.array(split_data.split(' '))
time = np.append(time, value[::2].astype(np.int))
data = np.append(data, value[1::2].astype(np.int))
# print(len(time), len(data))
xf = rfftfreq(N, 1 / SAMPLE_RATE)
points_per_freq = len(xf) / (SAMPLE_RATE / 2)
target_idx = int(points_per_freq * cut_off_freq)
# print(xf[target_idx])
yf = rfft(data)
yf[:target_idx] = 0
new_sig = np.round(irfft(yf, len(data)),3)
# print(time[-1])
# print(len(datas))
# print(len(yf),yf)
if threshold[0] == 'below':
th = new_sig < threshold[1]
elif threshold[0] == 'over':
th = new_sig > threshold[1]
elif threshold[0] == 'auto':
th = new_sig
threshold_edges = np.convolve([1, -1], th, mode='same')
thresholded_edge_indices = np.where(threshold_edges==1)[0]
filter_array = []
for idx, point in enumerate(thresholded_edge_indices):
if idx > 0:
if time[thresholded_edge_indices[idx]] - int(waveform[2]) >= time[thresholded_edge_indices[idx-1]]:
filter_array.append(True)
elif time[thresholded_edge_indices[idx]] - int(waveform[1]) >= time[thresholded_edge_indices[idx-1]]:
filter_array.append(True)
else:
filter_array.append(False)
else:
if time[thresholded_edge_indices[idx]] - int(waveform[1]) < 0:
filter_array.append(False)
else:
filter_array.append(True)
timemark_list = thresholded_edge_indices[filter_array]
if ret_data_type == 'all':
print(' '.join(str(x) for x in time))
print(' '.join(str(x) for x in new_sig))
print(' '.join(str(x) for x in timemark_list))
elif ret_data_type == 'partial':
return_sub_signal_lists = {}
for timemark in timemark_list:
start_index = timemark - floor(int(waveform[1] / DELTA_TIME))
end_index = timemark + ceil(int(waveform[2] / DELTA_TIME))
subSignal = new_sig[start_index: end_index + 2]
return_sub_signal_lists[str(timemark)] = ' '.join(str(x) for x in subSignal)
print(dumps(return_sub_signal_lists))
conn.commit()
conn.close()
exit