diff --git a/scripts/filters.py b/scripts/filters.py index a9daf50..763fca7 100644 --- a/scripts/filters.py +++ b/scripts/filters.py @@ -1,20 +1,35 @@ +import numpy as np from scipy.signal import butter, lfilter, decimate def anti_alias(data, sample_rate, max_frequency): + FILTER_HEADROOM = 1.2 nyquist_rate = 2 * max_frequency + downsample_factor = 1 if sample_rate > nyquist_rate: - filtered_data = lowpass_filter(data, max_frequency, sample_rate) + filtered_data = lowpass_filter(data, max_frequency * FILTER_HEADROOM, sample_rate) downsample_factor = int(sample_rate / nyquist_rate) - filtered_data = decimate(filtered_data, downsample_factor) sample_rate = sample_rate // downsample_factor else: filtered_data = data - return filtered_data, sample_rate + return filtered_data, sample_rate, downsample_factor + +def smoothing_filter(data, window_size=256): + window = np.ones(window_size) / window_size + return np.convolve(data, window, mode='same') + +# Stolen from selcald +def note(freq, length, amp=1.0, rate=44100): + if freq == 0: + data = np.zeros(int(length * rate)) + else: + t = np.linspace(0, length, int(length * rate)) + data = np.sin(2 * np.pi * freq * t) * amp + return data # These originally came from https://scipy.github.io/old-wiki/pages/Cookbook/ButterworthBandpass, # but they've been copied around the internet so many times that ChatGPT now produces them verbatim. diff --git a/scripts/selcal-detect.py b/scripts/selcal-detect.py index 72078a0..bcfafc8 100755 --- a/scripts/selcal-detect.py +++ b/scripts/selcal-detect.py @@ -5,85 +5,31 @@ import sys import numpy as np from scipy import signal from scipy.io import wavfile -from scipy.signal import butter, lfilter - - -tones = {} -with open('tones.csv', newline='') as csvfile: - reader = csv.DictReader(csvfile) - for row in reader: - tones[row['designator']] = float(row['frequency']) - -''' -def freq_of_key(midi_key): - return 440.0 * (2 ** ((midi_key - 69)/12)) - -tones = {} -for c in range(65, 90): - tones[c] = freq_of_key(c) -''' - -# Shamelessly lifted from -# https://scipy.github.io/old-wiki/pages/Cookbook/ButterworthBandpass -def butter_bandpass(lowcut, highcut, fs, order=5): - nyq = 0.5 * fs - low = lowcut / nyq - high = highcut / nyq - b, a = butter(order, [low, high], btype='band') - return b, a - -def butter_bandpass_filter(data, lowcut, highcut, fs, order=5): - b, a = butter_bandpass(lowcut, highcut, fs, order=order) - y = lfilter(b, a, data) - return y - -# tone synthesis -def note(freq, cycles, amp=32767.0, rate=44100): - len = cycles * (1.0/rate) - t = np.linspace(0, len, int(len * rate)) - if freq == 0: - data = np.zeros(int(len * rate)) - else: - data = np.sin(2 * np.pi * freq * t) * amp - return data.astype(int) - -def decimate_from_sample_rate(sample_rate): - if sample_rate == 44100: - return 4 # rate = 11025, Fmax = 5512.5 Hz - elif sample_rate == 48000: - return 5 # rate = 9600, Fmax = 4800 Hz - elif sample_rate == 22050: - return 2 # rate = 11025, Fmax = 5512.5 Hz - elif sample_rate == 11025: - return 1 # rate = 11025, Fmax = 5512.5 Hz - else: - raise ValueError("Sample rate not supported") +from tones import TONES +from filters import bandpass_filter, note, smoothing_filter, anti_alias +pure_sample_length = 0.1 if __name__ == '__main__': - # TODO JMT: What is this? - FLT_LEN = 2000 - file_name = sys.argv[1] sample_rate, data = wavfile.read(file_name) - print(f"{file_name}: {len(data)} samples @ {sample_rate} Hz") - decimate = decimate_from_sample_rate(sample_rate) - if decimate > 1: - data = signal.decimate(data, decimate) - sample_rate = sample_rate / decimate + if len(data.shape) == 2: + data = data.mean(axis=1) - print(f'Length after decimation: {len(data)} samples') + # Normalise + print(np.max(data)) + data = data / np.max(data) - data = butter_bandpass_filter(data, 270, 1700, sample_rate, order=8) + # TODO JMT: Find out why the correlation step fails when max frequency <= 2 * nyquist rate + data, sample_rate, decimation = anti_alias(data, sample_rate, 4800) + print(f'Length after decimation: {len(data)} samples (/{decimation}, {sample_rate})') - pure_signals = {tone:note(freq, FLT_LEN, rate=sample_rate) for tone,freq in tones.items()} + pure_signals = {tone:note(freq, pure_sample_length, rate=sample_rate) for tone,freq in TONES.items()} correlations = {tone:np.abs(signal.correlate(data, pure, mode='same')) for tone,pure in pure_signals.items()} + massaged = {tone:smoothing_filter(correlation) for tone,correlation in correlations.items()} - N = FLT_LEN // 8 # Rolling average length - cumsum_convolution = np.ones(N)/N - massaged = {tone:np.convolve(correlation, cumsum_convolution, mode='valid') for tone,correlation in correlations.items()} # Only import if we're actually plotting, these imports are pretty heavy. import pyqtgraph as pg @@ -100,7 +46,7 @@ if __name__ == '__main__': legend.setParentItem(legend_view) color_map = pg.colormap.get('CET-C6s') - colors = color_map.getLookupTable(nPts=len(tones)) + colors = color_map.getLookupTable(nPts=len(TONES)) for (tone, correlation), color in zip(massaged.items(), colors): line = plot.plot(correlation, pen=pg.mkPen(color=color), fillLevel=0.1, name=tone)