Here's what I have right now for fftfreq. I'm excited to have factored fftfreq out and added optional minimum and maximum frequency bounds. The rest of fourier.py is 1 email back. This interface does not facilitate the usecase of having minimum and maximum frequency bounds and simply desiring as many freq_count and freq_sample_rate as makes sense with the bounds. Thinking: the purpose of freq_sample_rate was to test against data that had been generated at sampled at a specific different rate, and was occurring in a larger recording. But maybe it would be more appropriate to consolidate that information into freq_count, min_freq, and max_freq, given in real recordings there is usually only one time base. # AGPL-3 Karl Semich 2022 import numpy as np def fftfreq(freq_count, sample_rate = None, min_freq = None, max_freq = None, dc_offset = True, freq_sample_rate = None): if sample_rate is None and freq_sample_rate is None: sample_rate = 1 freq_sample_rate = 1 elif sample_rate is None: sample_rate = freq_sample_rate elif freq_sample_rate is None: freq_sample_rate = sample_rate if not dc_offset: freq_count += 1 if min_freq is None: min_freq = freq_sample_rate / freq_count min_freq /= sample_rate if freq_count % 2 == 0: if max_freq is None: max_freq = freq_sample_rate / 2 max_freq /= sample_rate neg_freqs = np.linspace(-max_freq, -min_freq, num=freq_count//2, endpoint=True) pos_freqs = -neg_freqs[:0:-1] else: if max_freq is None: max_freq = freq_sample_rate * (freq_count - 1) / 2 / freq_count max_freq /= sample_rate pos_freqs = np.linspace(min_freq, max_freq, num=freq_count//2, endpoint=True) neg_freqs = -pos_freqs[::-1] return np.concatenate([ np.array([0] if dc_offset else []), pos_freqs, neg_freqs ]) def peak_pair_idcs(freq_data): freq_heights = abs(freq_data) # squares and sums the components paired_heights = freq_height[...,1:-1] + freq_height[...,2:] peak_idx = paired_heights.argmax(axis=-1, keepdims=True) + 1 return np.concatenate(peak_idx, peak_idx + 1, axis=-1)