Compute group delay of an audio file from STFT
The problem is that the stft function is splitting the signal up into different windows. That means that the signal from time $n$ to $n+N_{w}-1$ is multiplied by
$$ n, n+1, n+2, \ldots, n+N_w-1$$
instead of
$$
0, 1, 2, \ldots, N_w-1
$$
which is causing the scaling problem.

If I apply the group delay calculation from this derivation, I get:



where the top plot is the "data", the middle plot is the FFT of the data, and the bottom plot is the group delay. The group delay plot seems to be correct ($\tau_g = 49.5$, except for where the data is zero in the frequency domain.

To get this to work with the stft function, you're going to have to modulate the data with a sawtooth function rather than just np.arange(len(data)).

If I slice the sound file bachfugue.wav using the new, new code below I get the following image, which is more like what I'd expect.





New New Code Below

#!/opt/local/bin/python2.7
from scipy.signal import firwin
import numpy as np
import matplotlib.pyplot as plt
from scipy.io import wavfile


def group_delay(sig):
    b = np.fft.fft(sig)
    n_sig = np.multiply(sig, np.arange(len(sig)))
    br = np.fft.fft(n_sig)
    return np.divide(br, b + 0.01).real


rate, data = wavfile.read('bachfugue.wav')

N_w = 1024
L = np.floor(len(data)/N_w).astype(int)

full_delay = np.zeros((L, N_w))
for index in range(L):
    window = data[index:(index+N_w), 1].astype(float)
    full_delay[index][:] = group_delay(window)

plt.imshow(full_delay[0::10, 0::10], cmap=plt.get_cmap('hot'),aspect='auto')
plt.show()


New Code Below

#!/opt/local/bin/python2.7
from scipy.signal import firwin
import numpy as np
import matplotlib.pyplot as plt



N = 100
data = firwin(N, [0.05, 0.95], width=0.05, pass_zero=False)
n_data = np.multiply(data, np.arange(len(data)))
rate = 1

B = np.fft.fft(data)
Br = np.fft.fft(n_data)
group_delay2 = np.divide(Br, B + 0.01).real

plt.figure(1)
plt.subplot(311)
plt.plot(data)
plt.subplot(312)
plt.plot(np.abs(B))
plt.subplot(313)
plt.plot(group_delay2)

plt.show()


--

Original Answer

I've implemented it as I think it should be. I don't have your data, so I've just made a sinusoid.

The resulting plot is below. Code at end.

So I'm wondering what you should expect to see... more thought needed. I'll read that paper in more detail and see what I can make of it.





Now that I've had a chance to read through it, I'm wondering if that is correct.  Based on this derivation... I think the real way to do it might be different.  I'll try to code that up later.



Python Code

#!/opt/local/bin/python2.7
from scipy.signal import stft
import numpy as np
import matplotlib.pyplot as plt

N = 102400
rate = 1.0
data = []
for phi in ([2.0*3.14159265*x for x in range(0, N-1, 1)]):
    data.append(np.sin(phi))

n_data = np.multiply(data, np.arange(len(data)))
f, t, X = stft(data, rate, window="hamming", nperseg=800, return_onesided=False)
f_n, t_n, Y = stft(n_data, rate, window="hamming", nperseg=800, return_onesided=False)
# group_delay = (np.dot(X.real, Y.real) + np.dot(Y.imag, X.imag)) / np.power(np.abs(X), 2)

group_delay = np.divide(np.multiply(X.real, Y.real) + np.multiply(Y.imag, X.imag), np.power(np.abs(X), 2))

plt.imshow(group_delay)
plt.show()