raab/GP2023_chirp_detection
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

added more docstrings and moved chirpsim

Browse Source
This commit is contained in:
weygoldt 2023-01-12 08:50:28 +01:00
parent 21ff35c6ac
commit ce1c0fa319
4 changed files with 176 additions and 105 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

23
code/analysis.py
View File

@ -32,8 +32,8 @@ def main(folder):
for track_id in np.unique(ident):
# window_index for time array in time window
window_index = np.arange(len(idx))[(ident == track_id) &
(time[idx] >= t0) &
(time[idx] <= (t0+dt))]
(time[idx] >= t0) &
(time[idx] <= (t0+dt))]
freq_temp = freq[window_index]
time_temp = time[idx[window_index]]
#mean_freq = np.mean(freq_temp)
@ -46,20 +46,21 @@ def main(folder):
id = 10.
i = 10
window_index = np.arange(len(idx))[(ident == id) &
(time[idx] >= t0) &
(time[idx] <= (t0+dt))]
(time[idx] >= t0) &
(time[idx] <= (t0+dt))]
freq_temp = freq[window_index]
time_temp = time[idx[window_index]]
mean_freq = np.mean(freq_temp)
fdata = bandpass_filter(data_oi[:, i], rate=data.samplerate, lowf=mean_freq-5, highf=mean_freq+200)
fdata = bandpass_filter(
data_oi[:, i], rate=data.samplerate, lowf=mean_freq-5, highf=mean_freq+200)
fig, ax = plt.subplots()
ax.plot(np.arange(len(fdata))/data.samplerate, fdata, marker='*')
#plt.show()
#freqency analyis of filtered data
# plt.show()
# freqency analyis of filtered data
time_fdata = np.arange(len(fdata))/data.samplerate
roll_fdata = np.roll(fdata, shift=1)
period_index = np.arange(len(fdata))[(roll_fdata < 0) & (fdata>=0)]
period_index = np.arange(len(fdata))[(roll_fdata < 0) & (fdata >= 0)]
plt.plot(time_fdata, fdata)
plt.scatter(time_fdata[period_index], fdata[period_index], c='r')
@ -82,13 +83,10 @@ def main(folder):
# calculate the frequency
inst_freq = 1 / np.diff(true_zero)
filtered_inst_freq = gaussian_filter1d(inst_freq, 0.005)
fig, ax = plt.subplots()
fig, ax = plt.subplots()
ax.plot(filtered_inst_freq, marker='.')
# in 5 sekunden welcher fisch auf einer elektrode am
embed()
exit()
@ -98,7 +96,6 @@ def main(folder):
# fig, ax = plt.subplots(figsize=(20/2.54, 12/2.54))
# ax.plot(np.arange(len(data_oi[:, i])), data_oi[:, i])
pass

26
code/chirpdetection.py
View File

@ -121,7 +121,7 @@ def double_bandpass(
) -> tuple[np.ndarray, np.ndarray]:
"""
Apply a bandpass filter to the baseline of a signal and a second bandpass
filter above or below the baseline.
filter above or below the baseline, as specified by the search frequency.
Parameters
----------
@ -171,8 +171,13 @@ def main(datapath: str) -> None:
ident = np.load(datapath + "ident_v.npy", allow_pickle=True)
# set time window # <------------------------ Iterate through windows here
window_duration = 60 * data.samplerate
window_overlap = 0.3
window_duration = 60 * 5 # 5 minutes window
window_overlap = 30 # 30 seconds overlap
window_starts = np.arange(
time[0], time[-1], window_duration)
embed()
exit()
t0 = 3 * 60 * 60 + 6 * 60 + 43.5
dt = 60
@ -187,6 +192,14 @@ def main(datapath: str) -> None:
# <------------------------------------------ Find best electrodes here
# <------------------------------------------ Iterate through electrodes
# get indices for time array in time window
window_index = np.arange(len(idx))[
(ident == track_id) & (time[idx] >= t0) & (time[idx] <= (t0 + dt))
]
# filter baseline and above
freq_temp = freq[window_index]
time_temp = time[idx[window_index]]
electrode = 10
@ -224,14 +237,7 @@ def main(datapath: str) -> None:
# frequency where second filter filters
search_freq = -50
# get indices for time array in time window
window_index = np.arange(len(idx))[
(ident == track_id) & (time[idx] >= t0) & (time[idx] <= (t0 + dt))
]
# filter baseline and above
freq_temp = freq[window_index]
time_temp = time[idx[window_index]]
baseline, search = double_bandpass(
data_oi[:, electrode], data.samplerate, freq_temp, search_freq
)

153
code/modules/filters.py
View File

@ -2,96 +2,103 @@ from scipy.signal import butter, sosfiltfilt
import numpy as np
def bandpass_filter(data, rate, lowf=100, highf=1100):
def bandpass_filter(
data: np.ndarray,
rate: float,
lowf: float,
highf: float,
) -> np.ndarray:
"""Bandpass filter a signal.
Parameters
----------
data : np.ndarray
The data to be filtered
rate : float
The sampling rate
lowf : float
The low frequency cutoff
highf : float
The high frequency cutoff
Returns
-------
np.ndarray : The filtered data
"""
sos = butter(2, (lowf, highf), "bandpass", fs=rate, output="sos")
fdata = sosfiltfilt(sos, data)
return fdata
def highpass_filter(
data,
rate,
cutoff=100,
):
data: np.ndarray,
rate: float,
cutoff: float,
) -> np.ndarray:
"""Highpass filter a signal.
Parameters
----------
data : np.ndarray
The data to be filtered
rate : float
The sampling rate
cutoff : float
The cutoff frequency
Returns
-------
np.ndarray
The filtered data
"""
sos = butter(2, cutoff, "highpass", fs=rate, output="sos")
fdata = sosfiltfilt(sos, data)
return fdata
def lowpass_filter(
data,
rate,
cutoff=100,
):
data: np.ndarray,
rate: float,
cutoff: float
) -> np.ndarray:
"""Lowpass filter a signal.
Parameters
----------
data : np.ndarray
The data to be filtered
rate : float
The sampling rate
cutoff : float
The cutoff frequency
Returns
-------
np.ndarray
The filtered data
"""
sos = butter(2, cutoff, "lowpass", fs=rate, output="sos")
fdata = sosfiltfilt(sos, data)
return fdata
def envelope(data, rate, freq=100):
sos = butter(2, freq, "lowpass", fs=rate, output="sos")
envelope = np.sqrt(2) * sosfiltfilt(sos, np.abs(data))
return envelope
def envelope(data: np.ndarray, rate: float, freq: float) -> np.ndarray:
"""Calculate the envelope of a signal using a lowpass filter.
Parameters
----------
data : np.ndarray
The signal to calculate the envelope of
rate : float
The sampling rate of the signal
freq : float
The cutoff frequency of the lowpass filter
def create_chirp(
eodf=500,
chirpsize=100,
chirpduration=0.015,
ampl_reduction=0.05,
chirptimes=[0.05, 0.2],
kurtosis=1.0,
duration=1.0,
dt=0.00001,
):
"""create a fake fish eod that contains chirps at the given times. EOF is a simple sinewave. Chirps are modeled with Gaussian profiles in amplitude reduction and frequency ecxcursion.
Args:
eodf (int, optional): The chriping fish's EOD frequency. Defaults to 500 Hz.
chirpsize (int, optional): the size of the chrip's frequency excursion. Defaults to 100 Hz.
chirpwidth (float, optional): the duration of the chirp. Defaults to 0.015 s.
ampl_reduction (float, optional): Amount of amplitude reduction during the chrips. Defaults to 0.05, i.e. 5\%
chirptimes (list, optional): Times of chirp centers. Defaults to [0.05, 0.2].
kurtosis (float, optional): The kurtosis of the Gaussian profiles. Defaults to 1.0
dt (float, optional): the stepsize of the simulation. Defaults to 0.00001 s.
Returns:
np.ndarray: the time
np.ndarray: the eod
np.ndarray: the amplitude profile
np.adarray: tha frequency profile
Returns
-------
np.ndarray
The envelope of the signal
"""
p = 0.0
time = np.arange(0.0, duration, dt)
signal = np.zeros_like(time)
ampl = np.ones_like(time)
freq = np.ones_like(time)
ck = 0
csig = 0.5 * chirpduration / np.power(2.0 * np.log(10.0), 0.5 / kurtosis)
#csig = csig*-1
for k, t in enumerate(time):
a = 1.0
f = eodf
if ck < len(chirptimes):
if np.abs(t - chirptimes[ck]) < 2.0 * chirpduration:
x = t - chirptimes[ck]
gg = np.exp(-0.5 * np.power((x / csig) ** 2, kurtosis))
cc = chirpsize * gg
# g = np.exp( -0.5 * (x/csig)**2 )
f = chirpsize * gg + eodf
a *= 1.0 - ampl_reduction * gg
elif t > chirptimes[ck] + 2.0 * chirpduration:
ck += 1
freq[k] = f
ampl[k] = a
p += f * dt
signal[k] = -1 * a * np.sin(2 * np.pi * p)
return time, signal, ampl, freq
sos = butter(2, freq, "lowpass", fs=rate, output="sos")
envelope = np.sqrt(2) * sosfiltfilt(sos, np.abs(data))
return envelope

61
code/modules/simulations.py Normal file
View File

@ -0,0 +1,61 @@
import numpy as np
def create_chirp(
eodf=500,
chirpsize=100,
chirpduration=0.015,
ampl_reduction=0.05,
chirptimes=[0.05, 0.2],
kurtosis=1.0,
duration=1.0,
dt=0.00001,
):
"""create a fake fish eod that contains chirps at the given times. EOF is a simple sinewave. Chirps are modeled with Gaussian profiles in amplitude reduction and frequency ecxcursion.
Args:
eodf (int, optional): The chriping fish's EOD frequency. Defaults to 500 Hz.
chirpsize (int, optional): the size of the chrip's frequency excursion. Defaults to 100 Hz.
chirpwidth (float, optional): the duration of the chirp. Defaults to 0.015 s.
ampl_reduction (float, optional): Amount of amplitude reduction during the chrips. Defaults to 0.05, i.e. 5\%
chirptimes (list, optional): Times of chirp centers. Defaults to [0.05, 0.2].
kurtosis (float, optional): The kurtosis of the Gaussian profiles. Defaults to 1.0
dt (float, optional): the stepsize of the simulation. Defaults to 0.00001 s.
Returns:
np.ndarray: the time
np.ndarray: the eod
np.ndarray: the amplitude profile
np.adarray: tha frequency profile
"""
p = 0.0
time = np.arange(0.0, duration, dt)
signal = np.zeros_like(time)
ampl = np.ones_like(time)
freq = np.ones_like(time)
ck = 0
csig = 0.5 * chirpduration / np.power(2.0 * np.log(10.0), 0.5 / kurtosis)
#csig = csig*-1
for k, t in enumerate(time):
a = 1.0
f = eodf
if ck < len(chirptimes):
if np.abs(t - chirptimes[ck]) < 2.0 * chirpduration:
x = t - chirptimes[ck]
gg = np.exp(-0.5 * np.power((x / csig) ** 2, kurtosis))
cc = chirpsize * gg
# g = np.exp( -0.5 * (x/csig)**2 )
f = chirpsize * gg + eodf
a *= 1.0 - ampl_reduction * gg
elif t > chirptimes[ck] + 2.0 * chirpduration:
ck += 1
freq[k] = f
ampl[k] = a
p += f * dt
signal[k] = -1 * a * np.sin(2 * np.pi * p)
return time, signal, ampl, freq