Skip to content

Commit

Permalink
Merge pull request #17 from tillraab/master
Browse files Browse the repository at this point in the history 
changes an new script
  • Loading branch information
@jgrewe
jgrewe committed Mar 18, 2016
2 parents b4dd537 + 821f828 commit aa8ff3c
Show file tree
Hide file tree
Showing 6 changed files with 377 additions and 146 deletions.
55 changes: 45 additions & 10 deletions thunderfish/FishRecording.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -220,11 +220,21 @@ def best_window_algorithm(self, peak_no, mean_amplitudes, cvs, pks_th=0.15, ampl
valid_windows = valid_pks * valid_cv * valid_ampls

# If there is no best window, run the algorithm again with more flexible threshodlds.
if not True in valid_windows:
if not True in valid_windows and cvs_percentile_th == 100. and ampls_percentile_th == 0.:
print('\nWARNING. The recording %s seems to be of very bad quality for analysis. '
'Skipping recording...\n' % self._wavfile)
os.remove(self._wavfile)
quit()

elif not True in valid_windows:
print('\nNo best window found. Rerunning best_window_algorithm with more flexible arguments.\n')
if cvs_percentile_th <= 95.:
cvs_percentile_th += 5.
if ampls_percentile_th >= 5.:
ampls_percentile_th -= 5.
return self.best_window_algorithm(peak_no, mean_amplitudes, cvs,
ampls_percentile_th=ampls_percentile_th-5.,
cvs_percentile_th=cvs_percentile_th+5., plot_debug=plot_debug)
ampls_percentile_th=ampls_percentile_th,
cvs_percentile_th=cvs_percentile_th, plot_debug=plot_debug)
# This return is a Recursion! Need to return the value in the embeded function, otherwise the root_function
# will not return anything!

Expand Down Expand Up @@ -313,19 +323,44 @@ def plot_peaks_troughs(self, mod_filename):

pass

def type_detector(self, thres=.5):
def type_detector(self, thres=0.1, create_dataset=False, category = 'wave'):

pk_t, tr_t, _, _ = self.w_pt
pk_2_pk = pk_t[1:] - pk_t[:-1]
pk_2_pk = np.diff(pk_t)
pk_2_tr = np.abs(pk_t - tr_t)
med = np.median(2*pk_2_tr)

prop_in_2med = sum((pk_2_pk < 2*med) & (pk_2_pk > med))/float(len(pk_2_pk))
# in order to detect the type, we check the proportion of pk2pk time differences within 2* the median of pk2tr
# There should be a large proportion (~1.) for a wave type and a small proportion (~0.) for a pulse type.
min_n = np.min([len(pk_2_pk), len(pk_2_tr)])

r_tr = pk_2_tr[:min_n] / pk_2_pk[:min_n]
print np.min(r_tr), np.max(r_tr)
r_tr[r_tr>0.5] = 1 - r_tr[r_tr>0.5]
r_value = np.median(r_tr)
# r_value = np.median(pk_2_tr[:min_n] / pk_2_pk[:min_n])

if create_dataset:
if category is 'wave':
if not os.path.exists('wave_p2v_algor.npy'):
np.save('wave_p2v_algor.npy', np.array([]))
wave_p2v_algor = np.load('wave_p2v_algor.npy')
wave_p2v_algor = wave_p2v_algor.tolist()
wave_p2v_algor.append(r_value)
wave_p2v_algor = np.asarray(wave_p2v_algor)
np.save('wave_p2v_algor.npy', wave_p2v_algor)

elif category is 'pulse':
if not os.path.exists('pulse_p2v_algor.npy'):
np.save('pulse_p2v_algor.npy', np.array([]))
pulse_p2v_algor = np.load('pulse_p2v_algor.npy')
pulse_p2v_algor = pulse_p2v_algor.tolist()
pulse_p2v_algor.append(r_value)
pulse_p2v_algor = np.asarray(pulse_p2v_algor)
np.save('pulse_p2v_algor.npy', pulse_p2v_algor)

return 'pulse' if prop_in_2med < thres else 'wave'
else:
print("unknown fish category: %s" % category)
quit()

return 'pulse' if r_value < thres else 'wave', r_value
def plot_spectogram(self, ax):

fu_freq = self.fund_freq
Expand Down
123 changes: 80 additions & 43 deletions thunderfish/FishTracker.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@
import matplotlib.mlab as ml
from harmonic_tools import *
import sorting_tools as st
from IPython import embed


class FishTracker:
Expand Down Expand Up @@ -75,10 +76,13 @@ def processdata(self, data, fish_type, bwin, win_width, config_dict,
fishlists_dif_fresolution.append(fishlist)

fishlist = st.filter_fishes(fishlists_dif_fresolution)
# embed()
# filter fishlist so only 3 fishes with max strength are involved

psd_type = st.wave_or_pulse_psd(power_fres1, freqs_fres1,
if fish_type is 'pulse':
for i in np.arange(len(fishlist))[::-1]:
if fishlist[i][0][0] < 150:
fishlist.pop(i)

psd_type, mean_proportions = st.wave_or_pulse_psd(power_fres1, freqs_fres1,
data[(bwin * self.rate):(bwin * self.rate + win_width * self.rate)], self.rate,
self.fresolution)

Expand All @@ -95,7 +99,7 @@ def processdata(self, data, fish_type, bwin, win_width, config_dict,
# if psd_type is 'pulse' and fish_type is 'pulse':
# print 'HERE WE HAVE TO BUILD IN THE PULSEFISH ANALYSIS'
# # embed()
return power_fres1, freqs_fres1, psd_type, fish_type, fishlist
return power_fres1, freqs_fres1, psd_type, fish_type, fishlist, mean_proportions

def sort_my_wavefish(self):
"""
Expand Down Expand Up @@ -196,7 +200,7 @@ def exclude_short_files(self, data, index):
return good_file

def bw_psd_and_eod_plot(self, power, freqs, bwin, win_width, data, psd_type, fish_type, fishlist, pulse_data,
pulse_freq, output_path):
pulse_freq, output_path, mean_proportions, r_value, save_data_for_plot=False):
"""
Create figures showing the best window, its PSD and the the EOD of the fish
"""
Expand All @@ -222,12 +226,15 @@ def bw_psd_and_eod_plot(self, power, freqs, bwin, win_width, data, psd_type, fis
fs = 16 # fontsize of the axis labels
inch_factor = 2.54

fig, ax = plt.subplots(figsize=(plot_w / inch_factor, plot_h / inch_factor))
fig_all = plt.figure(figsize=(plot_w * 2 / inch_factor, plot_h * 2 / inch_factor))
ax1_all = fig_all.add_subplot(2, 2, (3, 4))
ax1_all = fig_all.add_subplot(2, 6, (7, 12))

power_dB = 10.0 * np.log10(power)
ax1_all.axis([0, 3000, min(power_dB[:len(freqs[freqs<3000])]), max(power_dB)+10])
ax1_all.plot(freqs, power_dB, lw=2, color='dodgerblue', alpha=0.7)

ax1_all.axis([0, 3000, -110, 0])
ax1_all.plot(freqs, 10.0 * np.log10(power), lw=2, color='dodgerblue', alpha=0.7)
if save_data_for_plot:
np.save('wave_psd_data.npy', np.array([freqs, power_dB]))

color = ['red', 'blue', 'green', 'cornflowerblue']
for color_no, fish in enumerate(ind):
Expand All @@ -245,45 +252,42 @@ def bw_psd_and_eod_plot(self, power, freqs, bwin, win_width, data, psd_type, fis
ax1_all.set_title('PSD of best window', fontsize=fs + 2)
sns.despine(fig=fig_all, ax=ax1_all, offset=10)

ax2_all = fig_all.add_subplot(2, 2, 2)
# ax2_all = fig_all.add_subplot(2, 6, (4, 6))

if psd_type is 'wave' or fish_type is 'wave':
fig3, ax3 = plt.subplots(figsize=(plot_w / inch_factor, plot_h / inch_factor))

#####################################################################################################
if fish_type is 'wave' and psd_type is 'wave':
ax2_all = fig_all.add_subplot(2, 6, (4, 6))
# fig3, ax3 = plt.subplots(figsize=(plot_w / inch_factor, plot_h / inch_factor))
# ToDo: Brasil data doesnt work properly: Line 260 Index out of bounds
eod_wtimes = np.arange(
len(data[(bwin * self.rate):(bwin * self.rate + round(self.rate * 1.0 / fishlist[ind[-1]][0][0] * 4))])
) * 1.0 / self.rate + bwin
eod_wampls = data[
(bwin * self.rate):(bwin * self.rate + round(self.rate * 1.0 / fishlist[ind[-1]][0][0] * 4))]
ax3.plot(eod_wtimes, eod_wampls, lw=2, color='dodgerblue', alpha=0.7,
label='Dominant frequency: %.2f Hz' % fishlist[ind[-1]][0][0])
ax3.tick_params(axis='both', which='major', labelsize=fs - 2)
plt.legend(loc='upper right')
ax3.set_xlabel('Time [sec]', fontsize=fs)
ax3.set_ylabel('Amplitude [a.u.]', fontsize=fs)
ax3.set_title('EOD-Waveform; %s' % self.filename, fontsize=fs + 2)
sns.despine(fig=fig3, ax=ax3, offset=10)
fig3.tight_layout()
fig3.savefig(plot_path + 'wave-EOD%.0f.pdf' % (len(glob.glob(plot_path + 'wave-EOD*.pdf')) + 1))
plt.close(fig3)

if fish_type is not 'pulse' or psd_type is not 'pulse':
ax2_all.plot(eod_wtimes, eod_wampls, lw=2, color='dodgerblue', alpha=0.7)
ax2_all.tick_params(axis='both', which='major', labelsize=fs - 2)
ax2_all.set_xlabel('Time [sec]', fontsize=fs)
ax2_all.set_ylabel('Amplitude [a.u.]', fontsize=fs)
ax2_all.set_title('EOD-Waveform', fontsize=fs + 2)
sns.despine(fig=fig_all, ax=ax2_all, offset=10)

# Soundtrace for a pulsefish-EOD #
# build mean and std over this data !
if psd_type is 'pulse' or fish_type is 'pulse':
# WHEN SINGLEWAVEFISH (len(fishlist) =1) gemittelter plot; anderenfall continue

ax2_all.plot(eod_wtimes, eod_wampls, lw=2, color='dodgerblue', alpha=0.7)
ax2_all.tick_params(axis='both', which='major', labelsize=fs - 2)
ax2_all.set_xlabel('Time [sec]', fontsize=fs)
ax2_all.set_ylabel('Amplitude [a.u.]', fontsize=fs)
ax2_all.set_title('EOD-Waveform', fontsize=fs + 2)
sns.despine(fig=fig_all, ax=ax2_all, offset=10)

if fish_type is 'pulse' or psd_type is 'pulse':
eod_plot_tw = 0.006
mean_pulse_data = []
std_pulse_data = []

for k in np.arange(len(pulse_data[1])):
datapoints_p_pulse = np.median([len(pulse_data[i]) for i in pulse_data.keys()])

for i in pulse_data.keys():
if len(pulse_data[i]) != int(datapoints_p_pulse):
# print i
del pulse_data[i]

for k in np.arange(int(datapoints_p_pulse)):

try:
tmp_mu = np.mean([pulse_data[pulse][k] for pulse in sorted(pulse_data.keys())])
Expand All @@ -301,8 +305,8 @@ def bw_psd_and_eod_plot(self, power, freqs, bwin, win_width, data, psd_type, fis

# get time for plot
plot_time = ((np.arange(len(mean_pulse_data)) * 1.0 / self.rate) - eod_plot_tw / 2) * 1000 # s to ms

if fish_type is 'pulse' and psd_type is 'pulse':
ax2_all = fig_all.add_subplot(2, 6, (4, 6))
ax2_all.plot(plot_time, mean_pulse_data, lw=2, color='dodgerblue', alpha=0.7, label='mean EOD')
ax2_all.plot(plot_time, up_std, lw=1, color='red', alpha=0.7, label='std EOD')
ax2_all.plot(plot_time, bottom_std, lw=1, color='red', alpha=0.7)
Expand All @@ -311,8 +315,41 @@ def bw_psd_and_eod_plot(self, power, freqs, bwin, win_width, data, psd_type, fis
ax2_all.set_ylabel('Amplitude [a.u.]', fontsize=fs)
ax2_all.set_title('Mean pulse-EOD', fontsize=fs + 2)
sns.despine(fig=fig_all, ax=ax2_all, offset=10)

text_ax = fig_all.add_subplot(2, 2, 1)
else:
ax3 = fig_all.add_subplot(2, 6, (3, 4))
ax3.plot(plot_time, mean_pulse_data, lw=2, color='dodgerblue', alpha=0.7, label='mean EOD')
ax3.plot(plot_time, up_std, lw=1, color='red', alpha=0.7, label='std EOD')
ax3.plot(plot_time, bottom_std, lw=1, color='red', alpha=0.7)
ax3.tick_params(axis='both', which='major', labelsize=fs - 2)
plt.legend(loc='upper right')
ax3.set_xlabel('Time [sec]', fontsize=fs)
ax3.set_ylabel('Amplitude [a.u.]', fontsize=fs)
ax3.set_title('EOD-Waveform pulse', fontsize=fs + 2)
sns.despine(fig=fig_all, ax=ax3, offset=10)
# fig3.tight_layout()
# fig3.savefig(plot_path + 'pulse-EOD%.0f.pdf' % (len(glob.glob(plot_path + 'pulse-EOD*.pdf')) + 1))
# plt.close(fig3)

ax5 = fig_all.add_subplot(2, 6, (5, 6))
eod_wtimes = np.arange(
len(data[(bwin * self.rate):(bwin * self.rate + round(self.rate * 1.0 / fishlist[ind[-1]][0][0] * 4))])
) * 1.0 / self.rate + bwin
eod_wampls = data[
(bwin * self.rate):(bwin * self.rate + round(self.rate * 1.0 / fishlist[ind[-1]][0][0] * 4))]

# WHEN SINGLEWAVEFISH (len(fishlist) =1) gemittelter plot; anderenfall continue

ax5.plot(eod_wtimes, eod_wampls, lw=2, color='dodgerblue', alpha=0.7)
ax5.tick_params(axis='both', which='major', labelsize=fs - 2)
ax5.set_xlabel('Time [sec]', fontsize=fs)
ax5.set_ylabel('Amplitude [a.u.]', fontsize=fs)
ax5.set_title('EOD-Waveform wave', fontsize=fs + 2)
sns.despine(fig=fig_all, ax=ax5, offset=10)

if psd_type is fish_type:
text_ax = fig_all.add_subplot(2, 6, (1, 3))
else:
text_ax = fig_all.add_subplot(2, 6, (1, 2))
text_ax.text(-0.1, 0.9, 'Filename:', fontsize=fs)
text_ax.text(0.5, 0.9, '%s' % self.filename, fontsize=fs)

Expand All @@ -329,10 +366,10 @@ def bw_psd_and_eod_plot(self, power, freqs, bwin, win_width, data, psd_type, fis
text_ax.text(0.5, 0.6, '%.2f Hz' % fishlist[ind[-1]][0][0], fontsize=fs)
else:
text_ax.text(0.5, 0.6, '%.2f Hz' % pulse_freq, fontsize=fs)
text_ax.text(0.5, 0.5, '%s' % fish_type, fontsize=fs)
text_ax.text(0.5, 0.5, '%s' % fish_type + ' (%.2f ; th = 0.1)' % r_value, fontsize=fs)

text_ax.text(-0.1, 0.4, 'PSD-Type:', fontsize=fs)
text_ax.text(0.5, 0.4, '%s' % psd_type, fontsize=fs)
text_ax.text(0.5, 0.4, '%s' % psd_type + ' (%.2f ; th = 0.27)' % mean_proportions, fontsize=fs)

text_ax.text(-0.1, 0.3, 'No. of wave-fishes:', fontsize=fs)
text_ax.text(0.5, 0.3, '%.0f' % len(fishlist), fontsize=fs)
Expand All @@ -359,9 +396,9 @@ def pulse_sorting(self, bwin, win_width, data):

# pulse frequency PROBLEMATISCH BEI MEHREREN PULSEFISHEN !!!
pulse_freq = len(th_time) / win_width
print ''
print 'Pulse-frequency:', pulse_freq
print ''
print('')
print('Pulse-frequency:', pulse_freq)
print('')

pulse_frequencies = 'fish_pulse.npy'
if not os.path.exists(pulse_frequencies):
Expand Down
2 changes: 1 addition & 1 deletion thunderfish/create_parallel_jobs.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,7 @@

# Create parallel jobs for user specific folder!
rec_dir = sys.argv[1]
job_name = 'parallel_jobs_for_' + sys.argv[-1].split('/')[-2] + '.txt'
job_name = 'parallel_jobs_for_' + sys.argv[-1].split('/')[-2] + '_' + sys.argv[-1].split('/')[-1] + '.txt'
if rec_dir[-1] != '/':
rec_dir += '/'

Expand Down
Loading

0 comments on commit aa8ff3c

Please sign in to comment.