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

adding peak detection

Browse Source
This commit is contained in:
wendtalexander 2023-01-11 15:49:42 +01:00
parent 98dee7a6ae
commit 9c2b4d0634
1 changed files with 29 additions and 5 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

32
code/chirpdetection.py
View File

@ -3,9 +3,12 @@ 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.stats import iqr
from scipy.signal import find_peaks
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 from modules.filters import bandpass_filter, envelope, highpass_filter
@ -32,7 +35,7 @@ def instantaneos_frequency(
true_zero = lower_time + lower_ratio * time_delta true_zero = lower_time + lower_ratio * time_delta
# create new time array # create new time array
inst_freq_time = true_zero[:-1] + 0.5 * np.diff(true_zero) inst_freq_time = true_zero[1:] + 0.5 * np.diff(true_zero)
# compute frequency # compute frequency
inst_freq = gaussian_filter1d(1 / np.diff(true_zero), 5) inst_freq = gaussian_filter1d(1 / np.diff(true_zero), 5)
@ -139,7 +142,7 @@ def main(datapath: str) -> None:
# track_id = ids # track_id = ids
# frequency where second filter filters # frequency where second filter filters
search_freq = 50 search_freq = -50
# get indices for time array in time window # get indices for time array in time window
window_index = np.arange(len(idx))[ window_index = np.arange(len(idx))[
@ -198,14 +201,35 @@ def main(datapath: str) -> None:
inst_freq_filtered = bandpass_filter( inst_freq_filtered = bandpass_filter(
baseline_freq, data.samplerate, lowf=15, highf=8000 baseline_freq, data.samplerate, lowf=15, highf=8000
) )
axs[6].plot(baseline_freq_time, np.abs(inst_freq_filtered))
# plot filtered and rectified envelope # plot filtered and rectified envelope
axs[4].plot(np.arange(len(baseline)) / axs[4].plot(np.arange(len(baseline)) /
data.samplerate, baseline_envelope) data.samplerate, baseline_envelope)
axs[5].plot(np.arange(len(baseline)) / axs[5].plot(np.arange(len(baseline)) /
data.samplerate, search_envelope) data.samplerate, search_envelope)
axs[6].plot(baseline_freq_time, np.abs(inst_freq_filtered))
# detect peaks baseline_enelope
prominence = iqr(baseline_envelope)
baseline_peaks, _ = find_peaks(baseline_envelope, prominence=prominence)
axs[4].scatter((np.arange(len(baseline)) /
data.samplerate)[baseline_peaks], baseline_envelope[baseline_peaks], c="red")
# detect peaks search_envelope
search_peaks, _ = find_peaks(search_envelope, height=0.0001)
axs[5].scatter((np.arange(len(baseline)) /
data.samplerate)[search_peaks], search_envelope[search_peaks], c="red")
# detect peaks inst_freq_filtered
inst_freq_peaks, _ = find_peaks(np.abs(inst_freq_filtered), height=2)
axs[6].scatter(baseline_freq_time[inst_freq_peaks], np.abs(inst_freq_filtered)[inst_freq_peaks], c="red")
#
axs[0].set_title("Spectrogram") axs[0].set_title("Spectrogram")
axs[1].set_title("Fitered baseline instanenous frequency") axs[1].set_title("Fitered baseline instanenous frequency")
axs[2].set_title("Fitered baseline") axs[2].set_title("Fitered baseline")