No więc tak. Do sygnału z fotodiody zastosowałem filtr dolnoprzepustowy, tj.:
# File name "Filters.py"
import scipy.signal as ss
def filt(sig, sf, cf, btype='higphass'):
"""
:param sig: signal.
:param sf: sampling frequency.
:param cf: cut frequencies - array.
:param btype: bandpass type.
:return: bandpassed signal.
"""
if btype == 'higphass' or btype == 'lowpass':
b, a = ss.butter(3, Wn=cf/(0.5*sf), btype=btype, analog=0, output='ba')
return ss.filtfilt(b, a, sig)
elif btype == 'bandstop' or btype == 'bandpass':
b, a = ss.butter(3, Wn=(cf[0]/(0.5*sf), cf[1]/(0.5*sf)), btype=btype, analog=0, output='ba')
return ss.filtfilt(b, a, sig)
I jego użycie dla 40 Hz:
import IBD.ElectricalStimulation.Filters as filt
filtered = filt.filt(signal, 1024, 40, btype='lowpass')
py.plot(time_scale, filtered)
Co dało mi:
Po czym wyliczyłem pochodną dla całego sygnału, dla kroku = 1, a całość podniosłem do kwadratu, żeby uwypuklić wartości.
# Derivate signal.
step = 1
accuracy_range = 9
derivative = np.diff(filtered, n=step)
derivative = np.append(derivative, np.zeros(step)) ** 2
derivative[derivative > accuracy_range] = np.max(filtered)
derivative[derivative < accuracy_range] = 0
py.plot(time_scale, derivative)
Otrzymałem w rezultacie coś takiego:
Prawie jestem w domu, ale niestety nie udaje mi się "uwzględnić" wszystkich "eventów", tj. wszystkich "istotnych" załamków.
UPDATE
Udało mi się jeszcze zwiększyć dokładność oznaczania, manipulując wartością obcinanej częstotliwości i zmienną accuracy_range
. Jednak tak, czy siak niektóre eventy pozostają nieoznaczone :/
# Bandpass signal.
import IBD.ElectricalStimulation.Filters as filt
filtered = filt.filt(signal, 1024, 33, btype='lowpass')
py.plot(time_scale, filtered)
# Derivate signal.
step = 1
accuracy_range = 1
derivative = np.diff(filtered, n=step)
derivative = np.append(derivative, np.zeros(step)) ** 2
derivative[derivative > accuracy_range] = np.max(filtered)
derivative[derivative < accuracy_range] = 0
py.plot(time_scale, derivative)
UPDATE 2
Wymyśliłem jeszcze coś takiego:
def walk_on_the_bitch(sig, t, interval=1000):
"""
:param sig: Pre-processed signal, ie. filtered.
:param t: Threshold.
:param interval: Interval between next value check (ms);
should be a predictive interval between each event.
"""
last_value = 0
interval_flag = True
interval_iterator = 0
markers = np.zeros(np.size(sig))
for i in np.arange(np.size(sig)):
absolute = np.abs(last_value - sig[i])
last_value = sig[i]
if interval_flag:
if absolute > t:
markers[i] = np.max(sig)
interval_flag = False
else:
if interval_iterator == interval:
interval_flag = True
interval_iterator = 0
else:
interval_iterator += 1
return markers
py.plot(time_scale, walk_on_the_bitch(filtered, 0.02))
Działa niby idealnie. Nie podoba mi się jednak ta drabinka ifów. Będę musiał jeszcze jednak nad tym pomyśleć. Dzięki chłopaki za pobudzenie kory :)