diff --git a/thunderfish/eodanalysis.py b/thunderfish/eodanalysis.py index 435f6761..b2754050 100644 --- a/thunderfish/eodanalysis.py +++ b/thunderfish/eodanalysis.py @@ -230,7 +230,8 @@ def analyze_wave(eod, freq, n_harm=6): return meod, props, spec_data -def analyze_pulse(eod, eod_times, thresh_fac=0.01, fresolution=1.0, percentile=1.0): +def analyze_pulse(eod, eod_times, thresh_fac=0.01, min_dist=50.0e-6, + fresolution=1.0, percentile=1.0): """ Analyze the EOD waveform of a pulse-type fish. @@ -243,6 +244,8 @@ def analyze_pulse(eod, eod_times, thresh_fac=0.01, fresolution=1.0, percentile=1 List of times of detected EOD peaks. thresh_fac: float Set the threshold for peak detection to the maximum pulse amplitude times this factor. + min_dist: float + Minimum distance between peak and troughs of the pulse. fresolution: float The frequency resolution of the power spectrum of the single pulse. percentile: float @@ -278,7 +281,7 @@ def analyze_pulse(eod, eod_times, thresh_fac=0.01, fresolution=1.0, percentile=1 and amplitude normalized to largest postive peak. power: 2-D array The power spectrum of a single pulse. First column are the frequencies, - second column the power in decibel scaled to the pulse rate. + second column the power. intervals: 1-D array List of inter-EOD intervals with extreme values removed. """ @@ -326,9 +329,12 @@ def analyze_pulse(eod, eod_times, thresh_fac=0.01, fresolution=1.0, percentile=1 # find smaller peaks: peak_idx, trough_idx = detect_peaks(meod[:,1], max_ampl*thresh_fac) peak_l = np.sort(np.concatenate((peak_idx, trough_idx))) - # remove mutliple peaks: - peak_list = [prange[np.argmax(np.abs(meod[prange,1]))] - for prange in np.split(peak_l, np.where(np.diff(meod[peak_l,1]>0.0)!=0)[0]+1)] + # remove mutliple peaks that do not oszillate around zero: + peak_list = np.array([prange[np.argmax(np.abs(meod[prange,1]))] + for prange in np.split(peak_l, np.where(np.diff(meod[peak_l,1]>0.0)!=0)[0]+1)]) + # remove multiple peaks that are too close and too small: + ridx = [(k, k+1) for k in np.where((np.diff(meod[peak_list,0]) < min_dist) & (meod[peak_list[:-1],1]<0.1))] + peak_list = np.delete(peak_list, ridx) # find P1: p1i = np.where(peak_list == max_idx)[0][0] offs = 0 if p1i <= 2 else p1i - 2 @@ -350,10 +356,9 @@ def analyze_pulse(eod, eod_times, thresh_fac=0.01, fresolution=1.0, percentile=1 data = np.zeros(nn) data[nn//2-n:nn//2+n] = meod[max_idx-n:max_idx+n,1] power, freqs = psd(data, samplerate, fresolution) - # store and scale: - dbpower = np.zeros((len(freqs), 2)) - dbpower[:,0] = freqs - dbpower[:,1] = decibel(power/period/period) + ppower = np.zeros((len(freqs), 2)) + ppower[:,0] = freqs + ppower[:,1] = power # store properties: props = {} @@ -375,7 +380,7 @@ def analyze_pulse(eod, eod_times, thresh_fac=0.01, fresolution=1.0, percentile=1 props['meaninterval'] = np.mean(intervals) props['stdinterval'] = np.std(intervals, ddof=1) - return meod, props, peaks, dbpower, intervals + return meod, props, peaks, ppower, intervals def eod_waveform_plot(eod_waveform, peaks, ax, unit=None, diff --git a/thunderfish/thunderfish.py b/thunderfish/thunderfish.py index d5d71d39..c43ce002 100644 --- a/thunderfish/thunderfish.py +++ b/thunderfish/thunderfish.py @@ -19,7 +19,7 @@ from .bestwindow import clip_amplitudes, best_window_indices from .harmonicgroups import add_psd_peak_detection_config, add_harmonic_groups_config from .checkpulse import check_pulse_width, add_check_pulse_width_config, check_pulse_width_args -from .powerspectrum import plot_decibel_psd, multi_resolution_psd +from .powerspectrum import decibel, plot_decibel_psd, multi_resolution_psd from .harmonicgroups import harmonic_groups, harmonic_groups_args, psd_peak_detection_args, fundamental_freqs, fundamental_freqs_and_power, colors_markers, plot_harmonic_groups from .consistentfishes import consistent_fishes from .eodanalysis import eod_waveform, analyze_wave, analyze_pulse, eod_waveform_plot @@ -125,6 +125,8 @@ def output_plot(base_name, pulse_fish, inter_eod_intervals, ############ # plot psd + if len(spec_data) > 0 and len(spec_data[0]) > 0.0: + ax3.plot(spec_data[0][:,0], decibel(5.0*eod_props[0]['EODf']**2.0*spec_data[0][:,1]), 'k', lw=1, alpha=0.2) if len(fishlist) > 0: colors, markers = colors_markers() plot_harmonic_groups(ax3, fishlist, max_freq=max_freq, max_groups=12, sort_by_freq=True, @@ -133,7 +135,6 @@ def output_plot(base_name, pulse_fish, inter_eod_intervals, frameon=False, bbox_to_anchor=(1.0, 1.1), loc='upper left', title='EOD frequencies') plot_decibel_psd(ax3, psd_data[0][1], psd_data[0][0], max_freq=max_freq, color='blue') - #ax3.plot(spec_data[0][:,0], spec_data[0][:,1], 'r', lw=2) ax3.set_title('Powerspectrum (%d detected wave-fish)' % len(fishlist), y=1.05) ########## @@ -165,12 +166,13 @@ def output_plot(base_name, pulse_fish, inter_eod_intervals, ################ - if not usedax5 and eod_props[0]['type'] == 'pulse': + if not usedax5 and len(eod_props) > 0 and eod_props[0]['type'] == 'pulse': usedax5 = True - smax = np.max(spec_data[0][:,1]) - ax5.plot(spec_data[0][:,0], spec_data[0][:,1] - smax, 'b', lw=3) - ax5.set_xscale('log') + db = decibel(spec_data[0][:,1]) + smax = np.nanmax(db) + ax5.plot(spec_data[0][:,0], db - smax, 'b', lw=3) ax5.set_xlim(1.0, 10000.0) + ax5.set_xscale('log') ax5.set_ylim(-60.0, 1.0) ax5.set_xlabel('Frequency [Hz]') ax5.set_ylabel('Power [dB]') @@ -244,7 +246,7 @@ def thunderfish(filename, channel=0, save_csvs=False, save_plot=False, cfg.add_section('Waveform selection:') cfg.add('maximumClippedFraction', 0.01, '', 'Take waveform of the fish with the highest power only if the fraction of clipped signals is below this value.') cfg.add('maximumFirstHarmonicAmplitude', 2.0, '', 'Skip waveform of wave-type fish if the ampltude of the first harmonic is higher than this factor times the amplitude of the fundamental.') - cfg.add('maximumSecondHarmonicAmplitude', 0.2, '', 'Skip waveform of wave-type fish if the ampltude of the second harmonic is higher than this factor times the amplitude of the fundamental. That is, the waveform appears to have twice the frequency than the fundamental.') + cfg.add('maximumSecondHarmonicAmplitude', 0.8, '', 'Skip waveform of wave-type fish if the ampltude of the second harmonic is higher than this factor times the amplitude of the fundamental. That is, the waveform appears to have twice the frequency than the fundamental.') cfg.add('maximumRMSError', 0.05, '', 'Skip waveform of wave-type fish if the root-mean-squared error relative to the peak-to-peak amplitude is larger than this number.') # load configuration from working directory and data directories: @@ -335,50 +337,23 @@ def thunderfish(filename, channel=0, save_csvs=False, save_plot=False, print(' none') fishlists.append(fishlist) # filter the different fishlists to get a fishlist with consistent fishes: - filtered_fishlist = consistent_fishes(fishlists, df_th=cfg.value('frequencyThreshold')) + fishlist = consistent_fishes(fishlists, df_th=cfg.value('frequencyThreshold')) if verbose > 0: - if len(filtered_fishlist) > 0: + if len(fishlist) > 0: print('fundamental frequencies consistent in all power spectra:') - print(' ' + ' '.join(['%.1f' % freq[0, 0] for freq in filtered_fishlist])) + print(' ' + ' '.join(['%.1f' % freq[0, 0] for freq in fishlist])) else: print('no fundamental frequencies are consistent in all power spectra') - # remove multiples of pulse fish frequency from fishlist: - if pulse_fish: - fishlist = [] - for fish in filtered_fishlist: - eodf = fish[0,0] - n = np.round(eodf*pulse_period) - pulse_eodf = n/pulse_period - if verbose > 0: - print('check wavefish at %.1f Hz versus %d-th harmonic of pulse frequency %.1f Hz at resolution of %.1f Hz' - % (eodf, n, pulse_eodf, 0.5*n*deltaf)) - # TODO: make the threshold a parameter! - if np.abs(eodf-pulse_eodf) > 0.5*n*deltaf: - fishlist.append(fish) - else: - if verbose > 0: - print('removed frequency %.1f Hz, because it is multiple of pulsefish %.1f' % (eodf, pulse_eodf)) - else: - fishlist = filtered_fishlist - - # save fish frequencies and amplitudes: - if save_csvs and found_bestwindow: - if len(fishlist) > 0: - # write csv file with main EODf and corresponding power in dB of detected fishes: - csv_matrix = fundamental_freqs_and_power(fishlist, cfg.value('powerNHarmonics')) - csv_name = os.path.join(output_folder, outfilename + '-wavefish-eodfs.csv') - header = ['fundamental frequency (Hz)', 'power (dB)'] - write_csv(csv_name, header, csv_matrix) - #if pulse_fish: - # TODO: write frequency amd parameter of pulse-fish to -pulsefish_eodfs.csv - - # analyse eod waveform: + # analysis results: eod_props = [] mean_eods = [] spec_data = [] peak_data = [] - inter_eod_intervals = [] + intervals = [] + power_thresh = [] + + # analyse eod waveform of pulse-fish: if pulse_fish: mean_eod_window = 0.002 mean_eod, eod_times = \ @@ -386,13 +361,26 @@ def thunderfish(filename, channel=0, save_csvs=False, save_plot=False, percentile=cfg.value('pulseWidthPercentile'), th_factor=cfg.value('pulseWidthThresholdFactor'), start=-mean_eod_window, stop=mean_eod_window) - mean_eod, props, peaks, power, inter_eod_intervals = analyze_pulse(mean_eod, eod_times, - fresolution=minfres) + mean_eod, props, peaks, power, intervals = analyze_pulse(mean_eod, eod_times, + fresolution=minfres) # TODO: add config parameter to analyze_pulse mean_eods.append(mean_eod) spec_data.append(power) peak_data.append(peaks) eod_props.append(props) + power_thresh = np.zeros(power.shape) + power_thresh[:,0] = power[:,0] + power_thresh[:,1] = 5.0*props['EODf']**2.0 * power[:,1] + + if len(power_thresh) > 0: + n = len(fishlist) + for k, fish in enumerate(reversed(fishlist)): + df = power_thresh[1,0] - power_thresh[0,0] + hfrac = float(np.sum(fish[:,1] < power_thresh[np.array(fish[:,0]//df, dtype=int),1]))/float(len(fish[:,1])) + if hfrac >= 0.5: + fishlist.pop(n-1-k) + if verbose > 0: + print('removed frequency %.1f Hz, because %.0f%% of the harmonics where below pulsefish threshold' % (fish[0,0], 100.0*hfrac)) # analyse EOD waveform of all wavefish: powers = np.array([np.sum(fish[:cfg.value('powerNHarmonics'), 1]) @@ -428,10 +416,11 @@ def thunderfish(filename, channel=0, save_csvs=False, save_plot=False, fishlist = [] eod_props = [] mean_eods = [] - inter_eod_intervals = [] + intervals = [] - # write eod waveforms: + # write results to files: if save_csvs and found_bestwindow: + # eod waveforms: for i, (mean_eod, sdata, pdata) in enumerate(zip(mean_eods, spec_data, peak_data)): header = ['time (s)', 'mean', 'std'] if mean_eod.shape[1] > 3: @@ -440,7 +429,7 @@ def thunderfish(filename, channel=0, save_csvs=False, save_plot=False, header, mean_eod) if len(sdata)>0: if sdata.shape[1] == 2: - header = ['frequency (Hz)', 'power (dB)'] + header = ['frequency (Hz)', 'power'] else: header = ['harmonics', 'frequency (Hz)', 'amplitude', 'phase'] if sdata.shape[1] > 4: @@ -451,10 +440,17 @@ def thunderfish(filename, channel=0, save_csvs=False, save_plot=False, header = ['P', 'time (ms)', 'amplitude', 'relaampl'] write_csv(os.path.join(output_folder, outfilename + '-peaks-%d.csv' % i), header, pdata) + # wavefish frequencies and amplitudes: + if len(fishlist) > 0: + # write csv file with main EODf and corresponding power in dB of detected fishes: + csv_matrix = fundamental_freqs_and_power(fishlist, cfg.value('powerNHarmonics')) + csv_name = os.path.join(output_folder, outfilename + '-wavefish-eodfs.csv') + header = ['fundamental frequency (Hz)', 'power (dB)'] + write_csv(csv_name, header, csv_matrix) if save_plot or not save_csvs: output_plot(outfilename, pulse_fish, - inter_eod_intervals, raw_data, samplerate, idx0, idx1, fishlist, + intervals, raw_data, samplerate, idx0, idx1, fishlist, mean_eods, eod_props, peak_data, spec_data, unit, psd_data, cfg.value('powerNHarmonics'), True, 3000.0, output_folder, save_plot=save_plot, show_plot=not save_csvs) diff --git a/thunderfish/version.py b/thunderfish/version.py index 74d7938f..c6f8db61 100644 --- a/thunderfish/version.py +++ b/thunderfish/version.py @@ -1 +1 @@ -__version__='1.2.0' # see http://semver.org/ +__version__='1.4.0' # see http://semver.org/