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

switched chirpdetection

Browse Source
This commit is contained in:
weygoldt 2023-01-11 14:11:38 +01:00
parent ed3791bc0d
commit ba116aaffe
1 changed files with 145 additions and 106 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

241
code/chirpdetection.py
View File

@ -3,181 +3,220 @@ import os
import numpy as np import numpy as np
from IPython import embed from IPython import embed
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
from scipy.ndimage import gaussian_filter1d
from thunderfish.dataloader import DataLoader from thunderfish.dataloader import DataLoader
from thunderfish.powerspectrum import spectrogram, decibel from thunderfish.powerspectrum import spectrogram, decibel
from scipy.ndimage import gaussian_filter1d
from modules.filters import bandpass_filter, envelope, highpass_filter, lowpass_filter from modules.filters import bandpass_filter, envelope, highpass_filter
def plot_spectogramm(ax, signal: np.ndarray, sampelrate: float) -> None: def instantaneos_frequency(
spec_power, spec_freqs, spec_times = spectrogram( signal: np.ndarray, samplerate: int
signal, ratetime=sampelrate, freq_resolution=50, overlap_frac=0.2 ) -> tuple[np.ndarray, np.ndarray]:
)
ax.pcolormesh(spec_times, spec_freqs, decibel(spec_power), vmin=-100, vmax=-50)
ax.set_ylim(500, 1200)
# calculate instantaneos frequency with zero crossings
roll_signal = np.roll(signal, shift=1)
time_signal = np.arange(len(signal)) / samplerate
period_index = np.arange(len(signal))[(roll_signal < 0) & (signal >= 0)]
def double_bandpass( upper_bound = np.abs(signal[period_index])
data: DataLoader, samplerate, freqs: np.ndarray, search_freq: float lower_bound = np.abs(signal[period_index - 1])
): upper_time = np.abs(time_signal[period_index])
lower_time = np.abs(time_signal[period_index - 1])
q25, q75 = np.percentile(freqs, [25, 75]) # create ratios
if q75 - q25 < 5: lower_ratio = lower_bound / (lower_bound + upper_bound)
baseline = np.median(freqs)
q25, q75 = baseline - 2.5, baseline + 2.5 # appy to time delta
# filter Baseline time_delta = upper_time - lower_time
filtered_baseline = bandpass_filter(data, samplerate, lowf=q25, highf=q75) true_zero = lower_time + lower_ratio * time_delta
# filter search area
filtered_searched_freq = bandpass_filter( # create new time array
data, samplerate, lowf=q25 + search_freq, highf=q75 + search_freq inst_freq_time = true_zero[:-1] + 0.5 * np.diff(true_zero)
# compute frequency
inst_freq = gaussian_filter1d(1 / np.diff(true_zero), 5)
return inst_freq_time, inst_freq
def plot_spectrogram(axis, signal: np.ndarray, samplerate: float) -> None:
# compute spectrogram
spec_power, spec_freqs, spec_times = spectrogram(
signal,
ratetime=samplerate,
freq_resolution=50,
overlap_frac=0.2,
) )
return (filtered_baseline, filtered_searched_freq) axis.pcolormesh(
spec_times,
spec_freqs,
decibel(spec_power),
)
axis.set_ylim(200, 1200)
def instantaneos_frequency(signal: np.ndarray, samplerate: int):
time_fdata = np.arange(len(signal)) / samplerate def double_bandpass(
roll_fdata = np.roll(signal, shift=1) data: DataLoader, samplerate: int, freqs: np.ndarray, search_freq: float
) -> tuple[np.ndarray, np.ndarray]:
period_index = np.arange(len(signal))[(roll_fdata < 0) & (signal >= 0)] # compute boundaries to filter baseline
q25, q75 = np.percentile(freqs, [25, 75])
upper_bound = np.abs(signal[period_index]) # check if percentile delta is too small
lower_bound = np.abs(signal[period_index - 1]) if q75 - q25 < 5:
median = np.median(freqs)
q25, q75 = median - 2.5, median + 2.5
upper_times = np.abs(time_fdata[period_index]) # filter baseline
lower_times = np.abs(time_fdata[period_index - 1]) filtered_baseline = bandpass_filter(data, samplerate, lowf=q25, highf=q75)
lower_ratio = lower_bound / (lower_bound + upper_bound) # filter search area
upper_ratio = upper_bound / (lower_bound + upper_bound) filtered_search_freq = bandpass_filter(
data, samplerate, lowf=q25 + search_freq, highf=q75 + search_freq
)
time_delta = upper_times - lower_times return (filtered_baseline, filtered_search_freq)
true_zero = lower_times + time_delta * lower_ratio
inst_freq = 1 / np.diff(true_zero)
filtered_inst_freq = gaussian_filter1d(inst_freq, 5)
# create new time axis
inst_freq_time = true_zero[:-1] + 0.5 * np.diff(true_zero)
return (inst_freq_time, filtered_inst_freq, true_zero) def main(datapath: str) -> None:
# load raw file
file = os.path.join(datapath, "traces-grid1.raw")
data = DataLoader(file, 60.0, 0, channel=-1)
def main(datapath: str): # load wavetracker files
# get the data
file = os.path.join(datapath, "traces-grid.raw")
data = DataLoader(datapath, 60.0, 0, channel=-1)
# load wavetracke files
time = np.load(datapath + "times.npy", allow_pickle=True) time = np.load(datapath + "times.npy", allow_pickle=True)
freq = np.load(datapath + "fund_v.npy", allow_pickle=True) freq = np.load(datapath + "fund_v.npy", allow_pickle=True)
ident = np.load(datapath + "ident_v.npy", allow_pickle=True)
idx = np.load(datapath + "idx_v.npy", allow_pickle=True) idx = np.load(datapath + "idx_v.npy", allow_pickle=True)
ident = np.load(datapath + "ident_v.npy", allow_pickle=True)
# make the right window for snipping # set time window # <------------------------ Iterate through windows here
t0 = 3 * 60 * 60 + 6 * 60 + 43.5 t0 = 3 * 60 * 60 + 6 * 60 + 43.5
dt = 60 dt = 60
start_index = t0 * data.samplerate start_index = t0 * data.samplerate
stop_index = (t0 + dt) * data.samplerate stop_index = (t0 + dt) * data.samplerate
# get the window with th data # load region of interest of raw data file
data_oi = data[start_index:stop_index, :] data_oi = data[start_index:stop_index, :]
# interate over the individuals # iterate through all fish
# track_id = np.unique(ident)[0] for track_id in np.unique(ident[~np.isnan(ident)])[:2]:
# <------------------------------------------ Find best electrodes here
# <------------------------------------------ Iterate through electrodes
# index of the electrode
electrode = 10 electrode = 10
for track_id in np.unique(ident[~np.isnan(ident)])[:2]:
# initialize plot
fig, axs = plt.subplots( fig, axs = plt.subplots(
7, 1, figsize=(20 / 2.54, 12 / 2.54), constrained_layout=True, sharex=True 7, 1, figsize=(20 / 2.54, 12 / 2.54), constrained_layout=True, sharex=True
) )
plot_spectogramm(axs[0], data_oi[:, electrode], data.samplerate) # plot spectrogram
plot_spectrogram(axs[0], data_oi[:, electrode], data.samplerate)
# plot wavetracker tracks to spectrogram
# for track_id in np.unique(ident): # <---------- Find freq gaps later
# here
# for track_id in np.unique(ident): # # get indices for time array in time window
# # window_index for time array in time window, fish data for time window # window_index = np.arange(len(idx))[
# window_index = np.arange(len(idx))[(ident == track_id) & # (ident == track_id) &
# (time[idx] >= t0) & # (time[idx] >= t0) &
# (time[idx] <= (t0+dt))] # (time[idx] <= (t0 + dt))
# ]
# freq_temp = freq[window_index] # freq_temp = freq[window_index]
# time_temp = time[idx[window_index]] # time_temp = time[idx[window_index]]
# axs[0].plot(time_temp - t0, freq_temp)
# axs[0].set_ylim(500, 1200) # axs[0].plot(time_temp-t0, freq_temp, lw=2)
# # define gap height # axs[0].set_ylim(500, 1000)
# # frequency for searching the chirp above the one fish
# track_id = ids
# frequency where second filter filters
search_freq = 50 search_freq = 50
# get indices for time array in time window
window_index = np.arange(len(idx))[ window_index = np.arange(len(idx))[
(ident == track_id) & (time[idx] >= t0) & (time[idx] <= (t0 + dt)) (ident == track_id) & (time[idx] >= t0) & (time[idx] <= (t0 + dt))
] ]
# filter baseline and above
freq_temp = freq[window_index] freq_temp = freq[window_index]
time_temp = time[idx[window_index]] time_temp = time[idx[window_index]]
brought_baseline = bandpass_filter(
data_oi[:, electrode],
data.samplerate,
lowf=np.mean(freq_temp) - 5,
highf=np.mean(freq_temp + 200),
)
baseline, search = double_bandpass( baseline, search = double_bandpass(
data_oi[:, electrode], data.samplerate, freq_temp, search_freq data_oi[:, electrode], data.samplerate, freq_temp, search_freq
) )
# calculate and plot the instantaneos freq # plot waveform of filtered signal
time_baseline_freq, basline_freq, ture_zeros = instantaneos_frequency( axs[2].plot(np.arange(len(baseline)) / data.samplerate, baseline)
# plot instatneous frequency
# broad_baseline = bandpass_filter(data_oi[:, electrode], data.samplerate, lowf=np.mean(
# freq_temp)-5, highf=np.mean(freq_temp)+200)
baseline_freq_time, baseline_freq = instantaneos_frequency(
baseline, data.samplerate baseline, data.samplerate
) )
inst_freq_filtered = bandpass_filter( axs[1].plot(baseline_freq_time, baseline_freq)
basline_freq, data.samplerate, lowf=1, highf=100
) # plot waveform of filtered search signal
axs[6].plot(time_baseline_freq, np.abs(inst_freq_filtered), marker=".") axs[3].plot(np.arange(len(baseline)) / data.samplerate, search)
axs[1].plot(time_baseline_freq, basline_freq, marker=".")
# compute envelopes
cutoff = 25 cutoff = 25
baseline_envelope = envelope(baseline, data.samplerate, cutoff) baseline_envelope = envelope(baseline, data.samplerate, cutoff)
axs[2].plot(ture_zeros, np.zeros_like(ture_zeros), marker=".", c="red")
axs[2].plot(np.arange(len(baseline)) / data.samplerate, baseline, c="blue")
axs[2].plot( axs[2].plot(
np.arange(len(baseline)) / data.samplerate, baseline_envelope, c="orange" np.arange(len(baseline)) / data.samplerate, baseline_envelope, c="orange"
) )
search_envelope = envelope(search, data.samplerate, cutoff) search_envelope = envelope(search, data.samplerate, cutoff)
axs[3].plot(np.arange(len(baseline)) / data.samplerate, search) axs[3].plot(
np.arange(len(baseline)) / data.samplerate, search_envelope, c="orange"
axs[3].plot(np.arange(len(baseline)) / data.samplerate, search_envelope) )
# filter and rectify envelopes # highpass filter envelopes
cutoff = 5 cutoff = 5
filtered_baseline_envelope = highpass_filter( baseline_envelope = highpass_filter(
baseline_envelope, data.samplerate, cutoff=cutoff baseline_envelope, data.samplerate, cutoff=cutoff
) )
filtered_searched_envelope = highpass_filter( # search_envelope = highpass_filter(
search_envelope, data.samplerate, cutoff=cutoff # search_envelope, data.samplerate, cutoff=cutoff)
)
# filter the envelopes bandpass # envelopes of filtered envelope of filtered baseline
baseline_envelope = envelope(
np.abs(baseline_envelope), data.samplerate, cutoff)
filtered_baseline_envelope = envelope( # search_envelope = bandpass_filter(
np.abs(filtered_baseline_envelope), data.samplerate, freq=5 # search_envelope, data.samplerate, lowf=lowf, highf=highf)
)
axs[4].plot( # bandpass filter the instantaneous
np.arange(len(baseline)) / data.samplerate, filtered_baseline_envelope inst_freq_filtered = bandpass_filter(
) baseline_freq, data.samplerate, lowf=15, highf=8000
axs[5].plot(
np.arange(len(baseline)) / data.samplerate, filtered_searched_envelope
) )
axs[6].plot(baseline_freq_time, np.abs(inst_freq_filtered))
axs[0].set_title("Spectogramm")
axs[1].set_title("Instantaneos Frequency") # plot filtered and rectified envelope
axs[2].set_title("Filtered Baseline") axs[4].plot(np.arange(len(baseline)) /
axs[3].set_title("Filtered Searched") data.samplerate, baseline_envelope)
axs[4].set_title("Filtered Baseline Envelope") axs[5].plot(np.arange(len(baseline)) /
axs[5].set_title("Filtered Searched Envelope") data.samplerate, search_envelope)
axs[0].set_title("Spectrogram")
axs[1].set_title("Fitered baseline instanenous frequency")
axs[2].set_title("Fitered baseline")
axs[3].set_title("Fitered above")
axs[4].set_title("Filtered envelope of baseline envelope")
axs[5].set_title("Search envelope")
axs[6].set_title("Filtered absolute instantaneous frequency")
plt.show() plt.show()
if __name__ == "__main__": if __name__ == "__main__":
datapath = "data/2022-06-02-10_00/" datapath = "../data/2022-06-02-10_00/"
main(datapath) main(datapath)