fixed python3 import pathes

setup.py

@@ -4,6 +4,11 @@
 setup(name='thunderfish',
       version='0.5.0',     # see http://semver.org/
       packages=find_packages(exclude=['contrib', 'doc', 'tests*']),
+      entry_points={
+        'console_scripts': [
+            'thunderfish = thunderfish.thunderfish:main',
+            'fishfinder = thunderfish.fishfinder:main',
+        ]},
       description='Algorithms and scripts for analyzing recordings of e-fish electric fields.',
       author='Jan Benda, Juan F. Sehuanes, Till Raab, Joerg Henninger, Jan Grewe, Fabian Sinz',
       requires=['numpy', 'matplotlib', 'audioio']

tests/test_harmonicgroups.py

@@ -0,0 +1,31 @@
+from nose.tools import assert_true, assert_equal, assert_almost_equal
+import numpy as np
+import thunderfish.fakefish as ff
+import thunderfish.powerspectrum as ps
+import thunderfish.harmonicgroups as hg
+
+
+def test_harmonic_groups():
+
+    # generate data:
+    samplerate = 44100.0
+    df = 0.5
+    eodfs = np.array([123.0, 321.0, 666.0, 668.0])
+    fish1 = ff.generate_wavefish(eodfs[0], samplerate, duration=8.0, noise_std=0.01,
+                                 amplitudes=[1.0, 0.5, 0.2, 0.1, 0.05])
+    fish2 = ff.generate_wavefish(eodfs[1], samplerate, duration=8.0, noise_std=0.01,
+                                 amplitudes=[1.0, 0.7, 0.2, 0.1])
+    fish3 = ff.generate_wavefish(eodfs[2], samplerate, duration=8.0, noise_std=0.01,
+                                 amplitudes=[10.0, 5.0, 1.0])
+    fish4 = ff.generate_wavefish(eodfs[3], samplerate, duration=8.0, noise_std=0.01,
+                                 amplitudes=[6.0, 3.0, 1.0])
+    data = fish1 + fish2 + fish3 + fish4
+
+    # analyse:
+    psd_data = ps.psd(data, samplerate, fresolution=df)
+    groups = hg.harmonic_groups(psd_data[1], psd_data[0])[0]
+    fundamentals = hg.fundamental_freqs(groups)
+
+    # check:
+    assert_true(np.all(np.abs(eodfs-fundamentals) < df),
+                'harmonic_groups() did not correctly detect all fundamental frequencies')

tests/test_peakdetection.py

@@ -1,4 +1,4 @@
-from nose.tools import assert_true, assert_equal, assert_almost_equal
+from nose.tools import assert_true, assert_equal, assert_almost_equal, assert_raises
 import numpy as np
 import thunderfish.peakdetection as pd
 
@@ -17,7 +17,7 @@ def test_detect_peaks():
     n = pt_indices[0]
     data[0:n] = 0.1 + 0.9 * np.arange(0.0, n) / n
     up = False
-    for i in xrange(0, len(pt_indices) - 1):
+    for i in range(len(pt_indices) - 1):
         n = pt_indices[i + 1] - pt_indices[i]
         if up:
             data[pt_indices[i]:pt_indices[i + 1]] = np.arange(0.0, n) / n
@@ -36,17 +36,11 @@ def test_detect_peaks():
     threshold = 0.5
     min_thresh = 0.3
 
-    peaks, troughs = pd.detect_peaks(data, 0.0)
-    assert_true(np.all(peaks == np.array([])) and np.all(troughs == np.array([])),
-                "detect_peaks(data, threshold) did not handle zero threshold")
+    assert_raises(ValueError, pd.detect_peaks, data, 0.0)
 
-    peaks, troughs = pd.detect_peaks(data, -1.0)
-    assert_true(np.all(peaks == np.array([])) and np.all(troughs == np.array([])),
-                "detect_peaks(data, threshold) did not handle negative threshold")
+    assert_raises(ValueError, pd.detect_peaks, data, -1.0)
 
-    peaks, troughs = pd.detect_peaks(data, threshold, time[:len(time) / 2])
-    assert_true(np.all(peaks == np.array([])) and np.all(troughs == np.array([])),
-                "detect_peaks(data, threshold) did not handle wrong time array")
+    assert_raises(IndexError, pd.detect_peaks, data, threshold, time[:len(time) / 2])
 
     peaks, troughs = pd.detect_peaks(data, threshold)
     assert_true(np.all(peaks == peak_indices),
@@ -100,7 +94,7 @@ def test_detect_dynamic_peaks():
     n = pt_indices[0]
     data[0:n] = 0.1 + 0.9 * np.arange(0.0, n) / n
     up = False
-    for i in xrange(0, len(pt_indices) - 1):
+    for i in range(len(pt_indices) - 1):
         n = pt_indices[i + 1] - pt_indices[i]
         if up:
             data[pt_indices[i]:pt_indices[i + 1]] = np.arange(0.0, n) / n
@@ -119,40 +113,26 @@ def test_detect_dynamic_peaks():
     threshold = 0.5
     min_thresh = 0.3
 
-    peaks, troughs = pd.detect_dynamic_peaks(data, 0.0, min_thresh, 0.5, time,
+    assert_raises(ValueError, pd.detect_dynamic_peaks, data, 0.0, min_thresh, 0.5, time,
                                              pd.accept_peak_size_threshold)
-    assert_true(np.all(peaks == np.array([])) and np.all(troughs == np.array([])),
-                "detect_dynamic_peaks(data, threshold) did not handle zero threshold")
 
-    peaks, troughs = pd.detect_dynamic_peaks(data, -1.0, min_thresh, 0.5, time,
+    assert_raises(ValueError, pd.detect_dynamic_peaks, data, -1.0, min_thresh, 0.5, time,
                                              pd.accept_peak_size_threshold)
-    assert_true(np.all(peaks == np.array([])) and np.all(troughs == np.array([])),
-                "detect_dynamic_peaks(data, threshold) did not handle negative threshold")
 
-    peaks, troughs = pd.detect_dynamic_peaks(data, threshold, 0.0, 0.5, time,
+    assert_raises(ValueError, pd.detect_dynamic_peaks, data, threshold, 0.0, 0.5, time,
                                              pd.accept_peak_size_threshold)
-    assert_true(np.all(peaks == np.array([])) and np.all(troughs == np.array([])),
-                "detect_dynamic_peaks(data, threshold) did not handle zero min_thresh")
 
-    peaks, troughs = pd.detect_dynamic_peaks(data, threshold, -1.0, 0.5, time,
+    assert_raises(ValueError, pd.detect_dynamic_peaks, data, threshold, -1.0, 0.5, time,
                                              pd.accept_peak_size_threshold)
-    assert_true(np.all(peaks == np.array([])) and np.all(troughs == np.array([])),
-                "detect_dynamic_peaks(data, threshold) did not handle negative min_thresh")
 
-    peaks, troughs = pd.detect_dynamic_peaks(data, threshold, min_thresh, 0.0, time,
+    assert_raises(ValueError, pd.detect_dynamic_peaks, data, threshold, min_thresh, 0.0, time,
                                              pd.accept_peak_size_threshold)
-    assert_true(np.all(peaks == np.array([])) and np.all(troughs == np.array([])),
-                "detect_dynamic_peaks(data, threshold) did not handle zero tau")
 
-    peaks, troughs = pd.detect_dynamic_peaks(data, threshold, min_thresh, -1.0, time,
+    assert_raises(ValueError, pd.detect_dynamic_peaks, data, threshold, min_thresh, -1.0, time,
                                              pd.accept_peak_size_threshold)
-    assert_true(np.all(peaks == np.array([])) and np.all(troughs == np.array([])),
-                "detect_dynamic_peaks(data, threshold) did not handle negative tau")
 
-    peaks, troughs = pd.detect_dynamic_peaks(data, threshold, min_thresh, 0.5, time[:len(time) / 2],
+    assert_raises(IndexError, pd.detect_dynamic_peaks, data, threshold, min_thresh, 0.5, time[:len(time) / 2],
                                              pd.accept_peak_size_threshold)
-    assert_true(np.all(peaks == np.array([])) and np.all(troughs == np.array([])),
-                "detect_dynamic_peaks(data, threshold) did not handle wrong time array")
 
     peaks, troughs = pd.detect_dynamic_peaks(data, threshold, min_thresh, 0.5, time,
                                              pd.accept_peak_size_threshold)

thunderfish/__init__.py

@@ -1,2 +1 @@
 __all__ = ['dataloader', 'configfile', 'peakdetection', 'bestwindow', 'powerspectrum', 'harmonicgroups', 'checkpulse', 'consistentfishes', 'eodanalysis', 'fakefish']
-from thunderfish import *

thunderfish/bestwindow.py

@@ -16,9 +16,8 @@
 plot_best_window(): visualization of the algorithm used in best_window_indices().
 """
 
-import warnings
 import numpy as np
-import peakdetection as pkd
+from .peakdetection import percentile_threshold, detect_peaks, trim_to_peak
 
 
 def clip_amplitudes(data, win_indices, min_fac=2.0, nbins=20,
@@ -221,18 +220,16 @@ def best_window_indices(data, samplerate, single=True, win_size=1., win_shift=0.
 
     # too little data:
     if len(data) / samplerate <= win_size:
-        warnings.warn('no best window found: not enough data')
-        return 0, 0
+        raise UserWarning('no best window found: not enough data')
 
     # threshold for peak detection:
-    threshold = pkd.percentile_threshold(data, samplerate, win_shift,
-                                         th_factor=th_factor, percentile=percentile)
+    threshold = percentile_threshold(data, samplerate, win_shift,
+                                     th_factor=th_factor, percentile=percentile)
 
     # detect large peaks and troughs:
-    peak_idx, trough_idx = pkd.detect_peaks(data, threshold)
+    peak_idx, trough_idx = detect_peaks(data, threshold)
     if len(peak_idx) == 0 or len(trough_idx) == 0:
-        warnings.warn('best_window(): no peaks or troughs detected')
-        return 0, 0
+        raise UserWarning('best_window(): no peaks or troughs detected')
 
     # compute cv of intervals, mean peak amplitude and its cv:
     invalid_cv = 1000.0
@@ -246,7 +243,7 @@ def best_window_indices(data, samplerate, single=True, win_size=1., win_shift=0.
         # indices of peaks and troughs inside analysis window:
         pinx = (peak_idx >= wtinx) & (peak_idx <= wtinx + win_size_indices)
         tinx = (trough_idx >= wtinx) & (trough_idx <= wtinx + win_size_indices)
-        p_idx, t_idx = pkd.trim_to_peak(peak_idx[pinx], trough_idx[tinx])
+        p_idx, t_idx = trim_to_peak(peak_idx[pinx], trough_idx[tinx])
         # interval statistics:
         ipis = np.diff(p_idx)
         itis = np.diff(t_idx)
@@ -275,14 +272,11 @@ def best_window_indices(data, samplerate, single=True, win_size=1., win_shift=0.
 
     # check:
     if len(mean_ampl[mean_ampl > 0.0]) <= 0:
-        warnings.warn('no finite amplitudes detected')
-        return 0, 0
+        raise UserWarning('no finite amplitudes detected')
     if len(cv_interv[cv_interv < invalid_cv]) <= 0:
-        warnings.warn('no valid interval cv detected')
-        return 0, 0
+        raise UserWarning('no valid interval cv detected')
     if len(cv_ampl[cv_ampl < invalid_cv]) <= 0:
-        warnings.warn('no valid amplitude cv detected')
-        return 0, 0
+        raise UserWarning('no valid amplitude cv detected')
 
     # cost function:
     cost = w_cv_interv * cv_interv + w_cv_ampl * cv_ampl - w_ampl * mean_ampl
@@ -292,8 +286,8 @@ def best_window_indices(data, samplerate, single=True, win_size=1., win_shift=0.
     valid_win_idx = np.nonzero(cost <= thresh)[0]
     cidx0 = valid_win_idx[0]  # start of current window
     cidx1 = cidx0 + 1  # end of current window
-    win_idx0 = cidx0  # start of largest window
-    win_idx1 = cidx1  # end of largest window
+    win_idx0 = cidx0   # start of largest window
+    win_idx1 = cidx1   # end of largest window
     i = 1
     while i < len(valid_win_idx):  # loop through all valid window positions
         if valid_win_idx[i] == valid_win_idx[i - 1] + 1:
@@ -496,10 +490,10 @@ def best_window_args(cfg):
         title = "test sines"
         data += 0.01 * np.random.randn(len(data))
     else:
-        import dataloader as dl
+        from .dataloader import load_data
 
         print("load %s ..." % sys.argv[1])
-        data, rate, unit = dl.load_data(sys.argv[1], 0)
+        data, rate, unit = load_data(sys.argv[1], 0)
         title = sys.argv[1]
 
     # determine clipping amplitudes:

thunderfish/checkpulse.py

@@ -6,9 +6,9 @@
 """
 
 import numpy as np
-import peakdetection as pkd
 import matplotlib.pyplot as plt
 from matplotlib.patches import Rectangle
+from .peakdetection import percentile_threshold, detect_peaks, trim_to_peak
 
 
 def check_pulse_width(data, samplerate, win_size=0.5, th_factor=0.8, percentile=0.1,
@@ -49,7 +49,7 @@ def ratio(peak_idx, trough_idx):
         :return: peak-ratio (float). The median of (peak-trough) / (peak-peak)
         :return: r_tr (array). The distribution of (peak-trough) / (peak-peak)
         """
-        peaks, troughs = pkd.trim_to_peak(peak_idx, trough_idx)
+        peaks, troughs = trim_to_peak(peak_idx, trough_idx)
 
         # get times of peaks and troughs, pk_times need to be floats!
         pk_times = peaks / float(samplerate)  # Actually there is no need to divide by samplerate.
@@ -70,11 +70,11 @@ def ratio(peak_idx, trough_idx):
         print('Analyzing Fish-Type...')
 
     # threshold for peak detection:
-    threshold = pkd.percentile_threshold(data, samplerate, win_size,
-                                         th_factor=th_factor, percentile=percentile)
+    threshold = percentile_threshold(data, samplerate, win_size,
+                                     th_factor=th_factor, percentile=percentile)
 
     # detect large peaks and troughs:
-    peak_idx, trough_idx = pkd.detect_peaks(data, threshold)
+    peak_idx, trough_idx = detect_peaks(data, threshold)
 
     pr_pvt, pvt_dist = ratio(peak_idx, trough_idx)
     pr_tvp, tvp_dist = ratio(trough_idx, peak_idx)
@@ -258,29 +258,29 @@ def plot_psd_proportion(freqs, power, proportions, percentiles, pulse_fish,
 if __name__ == "__main__":
     print("\nChecking checkpulse module ...\n")
     import sys
-    import bestwindow as bw
-    import powerspectrum as ps
-    import fakefish as ff
+    from .bestwindow import best_window
+    from .powerspectrum import multi_resolution_psd
+    from .fakefish import generate_monophasic_pulses, generate_biphasic_pulses, generate_triphasic_pulses, generate_alepto
 
     # generate data:
     rate = 44100.0
     if len(sys.argv) < 2:
-        data = ff.generate_biphasic_pulses(80.0, rate, 8.0)
+        data = generate_biphasic_pulses(80.0, rate, 8.0)
     elif sys.argv[1] == '-w':
-        data = ff.generate_alepto(600.0, rate, 8.0)
+        data = generate_alepto(600.0, rate, 8.0)
     elif sys.argv[1] == '-m':
-        data = ff.generate_monophasic_pulses(80.0, rate, 8.0)
+        data = generate_monophasic_pulses(80.0, rate, 8.0)
     elif sys.argv[1] == '-b':
-        data = ff.generate_biphasic_pulses(80.0, rate, 8.0)
+        data = generate_biphasic_pulses(80.0, rate, 8.0)
     elif sys.argv[1] == '-t':
-        data = ff.generate_triphasic_pulses(80.0, rate, 8.0)
+        data = generate_triphasic_pulses(80.0, rate, 8.0)
     else:  # load data given by the user
-        import dataloader as dl
+        from .dataloader import load_data
 
         file_path = sys.argv[1]
         print("loading %s ...\n" % file_path)
-        rawdata, rate, unit = dl.load_data(sys.argv[1], 0)
-        data, _ = bw.best_window(rawdata, rate)
+        rawdata, rate, unit = load_data(sys.argv[1], 0)
+        data, _ = best_window(rawdata, rate)
 
     # draw figure with subplots:
     fig1, ax1 = plt.subplots(nrows=3, ncols=1, figsize=(8., 12.))
@@ -291,7 +291,7 @@ def plot_psd_proportion(freqs, power, proportions, percentiles, pulse_fish,
     plt.tight_layout()
 
     fig2, ax2 = plt.subplots()
-    psd_data = ps.multi_resolution_psd(data, rate)
+    psd_data = multi_resolution_psd(data, rate)
     psd_type, proportions = check_pulse_psd(psd_data[0], psd_data[1], verbose=1,
                                             plot_data_func=plot_psd_proportion, ax=ax2, fs=12)
     plt.tight_layout()

thunderfish/chirp.py

@@ -8,9 +8,9 @@
 
 import numpy as np
 import matplotlib.pyplot as plt
-import harmonicgroups as hg
-import powerspectrum as ps
-import peakdetection as pkd
+from .harmonicgroups import harmonic_groups
+from .powerspectrum import spectrogram
+from .peakdetection import std_threshold, detect_peaks, trim_to_peak
 
 
 def true_chirp_power_drop(chirp_time_idx, power, power_window=100):
@@ -99,9 +99,9 @@ def     chirp_detection(spectrum, freqs, time, fishlist=None, fundamentals=None,
         power_diff = np.diff(power)
 
         # peakdetection in the power_diff to detect drops in power indicating chrips
-        threshold = pkd.std_threshold(power_diff)
-        peaks, troughs = pkd.detect_peaks(power_diff, threshold)
-        troughs, peaks = pkd.trim_to_peak(troughs, peaks) # reversed troughs and peaks in output and input to get trim_to_troughs
+        threshold = std_threshold(power_diff)
+        peaks, troughs = detect_peaks(power_diff, threshold)
+        troughs, peaks = trim_to_peak(troughs, peaks) # reversed troughs and peaks in output and input to get trim_to_troughs
 
         # exclude peaks and troughs with to much time diff to be a chirp
         # ToDO: not nice !!!
@@ -159,7 +159,7 @@ def chirp_detection_plot(enu, chirp_time, time, power, power2, power_diff, funda
     plt.legend(loc='upper right', bbox_to_anchor=(1, 1), frameon=False)
 
 
-def chirp_analysis(data, samplerate, cfg):
+def chirp_analysis(data, samplerate):
     """
     Performs all steps to detect chirps in a given dataset. This includes spectrogram calculation, fish detection and
     analysing of specific frequency bands.
@@ -168,14 +168,13 @@ def chirp_analysis(data, samplerate, cfg):
 
     :param data: (array) data.
     :param samplerate: (float) smaplerate of the data.
-    :param cfg:(dict) HAS TO BE REMOVED !!!!
     :param min_power: (float) minimal power of the fish fundamental to include this fish in chirp detection.
     """
-    spectrum, freqs, time = ps.spectrogram(data, samplerate, fresolution=2., overlap_frac=0.95)
+    spectrum, freqs, time = spectrogram(data, samplerate, fresolution=2., overlap_frac=0.95)
 
     power = np.mean(spectrum, axis=1) # spectrum[:, t0:t1] to only let spectrum of certain time....
 
-    fishlist = hg.harmonic_groups(freqs, power, cfg)[0]
+    fishlist = harmonic_groups(freqs, power)[0]
 
     chirp_time, chirp_freq = chirp_detection(spectrum, freqs, time, fishlist, plot_data_func=chirp_detection_plot)
 
@@ -190,14 +189,11 @@ def chirp_analysis(data, samplerate, cfg):
     # '2016_04_27__downstream_stonewall_at_pool' made in colombia, 2016.
     ###
     import sys
-    import dataloader as dl
-    import config_tools as ct
+    from .dataloader import load_data
 
-    cfg = ct.get_config_dict()
+    data_file = sys.argv[1]
+    raw_data, samplerate, unit = load_data(data_file, channel=0)
 
-    audio_file = sys.argv[1]
-    raw_data, samplerate, unit = dl.load_data(audio_file, channel=0)
+    chirp_time, chirp_freq = chirp_analysis(raw_data, samplerate)
 
-    chirp_time, chirp_freq = chirp_analysis(raw_data, samplerate, cfg)
-
-    # power = np.mean(spectrum[:, t:t + nffts_per_psd], axis=1)
+    # power = np.mean(spectrum[:, t:t + nffts_per_psd], axis=1)