-
Notifications
You must be signed in to change notification settings - Fork 9
/
Copy pathgroup_delay.py
100 lines (86 loc) · 3.26 KB
/
group_delay.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
# select time range to correlate
# select number of bands
# plot x: freq, y: phase delay, correlation
import logging
from util.correlation import xcorr, parabolic
logging.basicConfig(level=logging.INFO)
# logging.basicConfig(level=logging.DEBUG)
import numpy as np
import soundfile as sf
from dropouts_gui import pairwise
import matplotlib.pyplot as plt
from util import filters
def get_group_delay(ref_sig, src_sig):
# first get the time range for both
# apply bandpass
# split into pieces and look up the delay for each
# correlate all the pieces
# sr = self.sr
sr = 44100
f_upper = 2000
f_lower = 10
# num = 10
bandwidth = 45
num_bands = int((f_upper-f_lower) / bandwidth)
t0 = 153.2
t1 = 156.0
min_corr = 0.6
s_start = int(t0 * sr)
s_end = int(t1 * sr)
s_dur = s_end-s_start
s_window = np.hanning(s_dur)
# or logspace?
# band_limits = np.linspace(f_lower, f_upper, num_bands)
band_limits = np.logspace(np.log2(f_lower), np.log2(f_upper), num=num_bands, endpoint=True, base=2, dtype=np.uint16)
lags = []
correlations = []
band_centers = []
for f_lower_band, f_upper_band in pairwise(band_limits):
logging.info(f"lower {f_lower_band:.1f}, upper {f_upper_band:.1f}, width {f_upper_band-f_lower_band:.1f}")
ref_s = filters.butter_bandpass_filter(ref_sig[s_start:s_end], f_lower_band, f_upper_band, sr, order=1)
src_s = filters.butter_bandpass_filter(src_sig[s_start:s_end], f_lower_band, f_upper_band, sr, order=1)
# plt.plot(ref_s, label="ref_s")
# plt.plot(src_s, label="src_s")
# res = np.correlate(ref_s*s_window, src_s*s_window, mode="same")
# res = scipy.signal.correlate(ref_s*s_window, src_s*s_window, mode='same', method='auto')
# res = scipy.signal.correlate(ref_s, src_s, mode='same', method='auto')
# res = xcorr(ref_s*s_window, src_s*s_window, mode='same')
res = xcorr(ref_s, src_s, mode='same')
# this should maybe hanned before argmax to kill obvious outliers
i_peak = np.argmax(res*s_window)
# interpolate the most accurate fit
i_interp, corr = parabolic(res, i_peak)
v = (s_dur // 2) - i_interp
if corr > min_corr:
lags.append(v)
correlations.append(corr)
band_center = (f_lower_band+f_upper_band)/2
band_centers.append(band_center)
else:
logging.warning(f"Band had too little correlation {corr}")
#
# plt.title('Digital filter group delay')
# plt.plot(ref_s)
# plt.plot(src_s)
# plt.ylabel('a')
# plt.xlabel('t')
# plt.show()
plot_corr_lag(band_centers, correlations, lags)
def plot_corr_lag(band_centers, correlations, lags):
fig, ax1 = plt.subplots()
color = 'tab:red'
ax1.set_xlabel('frequency [Hz]')
ax1.set_ylabel('lag [sample]', color=color)
ax1.plot(band_centers, lags, color=color)
ax1.tick_params(axis='y', labelcolor=color)
ax2 = ax1.twinx() # instantiate a second axes that shares the same x-axis
color = 'tab:blue'
ax2.set_ylabel('correlation', color=color) # we already handled the x-label with ax1
ax2.plot(band_centers, correlations, color=color)
ax2.tick_params(axis='y', labelcolor=color)
fig.tight_layout() # otherwise the right y-label is slightly clipped
plt.show()
fp = "C:/Users/arnfi/Music/The Beatles/Revolver Companion/Rain/rain rhythm sinc az.wav"
ref_ob = sf.SoundFile(fp)
ref_sig = ref_ob.read(always_2d=True, dtype='float32')
get_group_delay(ref_sig[:, 0], ref_sig[:, 1])