信息处理 - 从 STFT 计算音频文件的群延迟 - 吾爱随笔录

从 STFT 计算音频文件的群延迟

信息处理过滤器声音的软件实现群延迟音频处理

2022-02-09 12:09:31

我想重现这篇论文的结果，其中作者使用群延迟图作为机器学习模型的输入。机器学习方面不是我遇到的问题。我在计算音频文件的群延迟函数时遇到问题。

首先，我知道群延迟本身是滤波器的一个属性。然而，该论文的作者写了这样的东西

群时延 [10] 定义为 STFT 相位谱的负导数：由于等式 (2) 的实现需要解缠相位谱，群延迟函数也可以仅使用幅度值来计算：其中和分别表示实部和虚部。和分别表示和的 STFT 。
$τ (ω, t) = - \frac{d (θ (ω, t))}{d ω} (2)$ $\tau(\omega,t)=-\frac{d(\theta(\omega,t))}{d\omega}\quad(2)$ $τ (ω, t) = \frac{X_{R} (ω, t) Y_{R} (ω, t) + Y_{I} (ω, t) X_{I} (ω, t)}{| X (ω, t) |^{2}} (3)$ $\tau(\omega,t)=\frac{X_{R}(\omega,t)Y_{R}(\omega,t)+Y_{I}(\omega,t)X_{I}(\omega,t)}{|X(\omega,t)|^2}\quad(3)$ $R$ $I$ $X(\omega,t)$ $Y(\omega,t)$ $x(n)$ $nx(n)$

我相信，这应该会在论文中生成如图 3 所示的图像，看起来像这样（在沿时间轴截断或填充到 256 的长度之后）

和都是非对称的，因此无法进行矩阵乘法。这是原始论文，在本文和 [10] 中都被引用，作者在其中推导方程 (3)。 $X$ $Y$

据我了解，矩阵乘法是标准乘法，而不是元素乘法（元素乘法不会产生论文中的结果）。在下面的实现中，我使用长度为 800 个样本的窗口，因为本文的作者使用 50 ms 窗口，这对应于 16 kHz 采样率音频文件的 800 个样本。我试图在python中实现它，如下所示：

from scipy.io import wavfile
from scipy.signal import stft, get_window
import numpy as np
import matplotlib.pyplot as plt

rate, data = wavfile.read("test.wav")

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)

但代码失败了

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-8-0058b9785caf> in <module>
----> 1 group_delay = (np.dot(X.real, Y.real) + np.dot(Y.imag, X.imag)) / np.power(np.abs(X), 2)

ValueError: shapes (800,109) and (800,109) not aligned: 109 (dim 1) != 800 (dim 0)

这是相当明显的。我也尝试在 Matlab 中实现它，希望这可能是 python 错误，但它也失败了。

[x,fs] = audioread('test.wav');
n_x = x .* [1:size(x,1)];
X = spectrogram(x, 800, 400, 2048, fs);
Y = spectrogram(n_x, 800, 400, 2048, fs);
T = (real(X) * real(Y) + imag(Y) * imag(X)) / abs(X)^2;

如何使用上述方程计算群延迟函数

1个回答

问题是该stft功能将信号分成不同的窗口。这意味着从时间到的信号乘以而不是导致缩放问题。 $n$ $n+N_{w}-1$

n, n + 1, n + 2, \dots, n + N_{w} - 1

$n, n+1, n+2, \ldots, n+N_w-1$

0, 1, 2, \dots, N_{w} - 1

$0, 1, 2, \ldots, N_w-1$

如果我从这个推导应用群延迟计算，我得到：

其中顶部图是“数据”，中间图是数据的 FFT，底部图是群延迟。群延迟图似乎是正确的（，但频域中的数据为零。 $\tau_g = 49.5$

为了让它与stft函数一起工作，你将不得不用锯齿函数而不是仅仅调制数据np.arange(len(data))。

如果我使用下面的新代码对声音文件进行切片bachfugue.wav，我会得到以下图像，这更像我所期望的。

下面的新代码

#!/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()

下面的新代码

#!/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()

原始答案

我已经按照我认为应该的方式实现了它。我没有你的数据，所以我只是做了一个正弦曲线。

结果图如下。代码在最后。

所以我想知道你应该期待看到什么......需要更多的思考。我会更详细地阅读那篇论文，看看我能做些什么。

现在我有机会通读它，我想知道这是否正确。基于这个推导......我认为真正的方法可能会有所不同。稍后我会尝试编写代码。

Python代码

#!/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()

其它你可能感兴趣的问题