文章代码仓库：https://gitee.com/LXP-Never/window_fun
窗函数贯穿整个语音信号处理，语音信号是一个非平稳的时变信号，但“**短时间内可以认为语音信号是平稳时不变的，一般 10~30ms**。
对连续的语音分帧做STFT处理，等价于截取一段时间信号，对其进行周期性延拓，从而变成无限长序列，并对该无限长序列做FFT变换，这一截断并不符合傅里叶变换的定义。因此，会导致频谱泄漏和混叠。
频谱泄漏：如果不加窗，默认就是矩形窗，时域的乘积就是频域的卷积，使得频谱以实际频率值为中心， 以窗函数频谱波形的形状向两侧扩散，指某一频点能量扩散到相邻频点的现象，会导致幅度较小的频点淹没在幅度较大的频点泄漏分量中
频谱混叠：会在分段拼接处引入虚假的峰值，进而不能获得准确的频谱情况
加窗的目的：让一帧信号的幅度在两端渐变到 0，渐变对傅里叶变换有好处，可以让频谱上的各个峰更细，不容易糊在一起，从而减轻频谱泄漏和混叠的影响。
加窗的代价：一帧信号两端的部分被削弱了，没有像中央的部分那样得到重视。弥补的办法就是相互重叠。相邻两帧的起始位置的时间差叫做帧移，常见的取法是取为帧长的一半。
对于语音，窗函数常选汉宁窗(Hanning)、汉明窗(Hamming)、sqrthann及其改进窗，他们的时域波形和幅频响应如下所示：
1、汉宁窗(Hann)
$$w(n) = 0.5 - 0.5 \cos\left(\frac{2\pi{n}}{M-1}\right) \qquad 0 \leq n \leq M-1$$
2、汉明窗(Hamming)
$$w(n) = 0.54 - 0.46 \cos\left(\frac{2\pi{n}}{M-1}\right) \qquad 0 \leq n \leq M-1$$
# -*- coding:utf-8 -*-
# Author:凌逆战 | Never
# Date: 2023/1/1
"""
绘制 窗函数和对应的频率响应
"""
import numpy as np
from numpy.fft import rfft
import matplotlib.pyplot as plt
window_len = 60
# frequency response
def frequency_response(window, window_len=window_len, NFFT=2048):
A = rfft(window, NFFT) / (window_len / 2) # (513,)
mag = np.abs(A)
freq = np.linspace(0, 0.5, len(A))
# 忽略警告
with np.errstate(divide='ignore', invalid='ignore'):
response = 20 * np.log10(mag)
response = np.clip(response, -150, 150)
return freq, response
def Rectangle_windows(win_length):
# 矩形窗
return np.ones((win_length))
def Voibis_windows(win_length):
""" Voibis_windows窗函数，RNNoise使用的是它，它满足Princen-Bradley准则。
:param x:
:param win_length: 窗长
:return:
"""
x = np.arange(0, win_length)
return np.sin((np.pi / 2) * np.sin((np.pi * x) / win_length) ** 2)
def sqrt_hanning_windows(win_length, mode="periodic"):
# symmetric: 对称窗，主要用于滤波器的设计
# periodic: 周期窗，常用于频谱分析
if mode == "symmetric":
haning_window = np.hanning(win_length)
sqrt_haning_window = np.sqrt(haning_window)
elif mode == "periodic":
haning_window = np.hanning(win_length+1)
sqrt_haning_window = np.sqrt(haning_window)
sqrt_haning_window = sqrt_haning_window[0:-1].astype('float32')
return sqrt_haning_window
Rectangle_windows = Rectangle_windows(window_len)
hanning_window = np.hanning(M=window_len)
print(np.argmax(hanning_window))
sqrt_hanning_windows = sqrt_hanning_windows(window_len)
hamming_window = np.hamming(M=window_len)
Voibis_windows = Voibis_windows(window_len)
blackman_window = np.blackman(M=window_len)
bartlett_window = np.bartlett(M=window_len)
kaiser_window = np.kaiser(M=window_len, beta=14)
plt.figure()
plt.plot(Rectangle_windows, label="Rectangle")
plt.plot(hanning_window, label="hanning")
plt.plot(sqrt_hanning_windows, label="sqrt_hanning")
plt.plot(hamming_window, label="hamming")
plt.plot(Voibis_windows, label="Voibis")
plt.plot(blackman_window, label="blackman")
plt.plot(bartlett_window, label="bartlett")
plt.plot(kaiser_window, label="kaiser")
plt.legend()
plt.tight_layout()
plt.show()
freq, Rectangle_FreqResp = frequency_response(Rectangle_windows, window_len)
freq, hanning_FreqResp = frequency_response(hanning_window, window_len)
freq, sqrt_hanning_FreqResp = frequency_response(sqrt_hanning_windows, window_len)
freq, hamming_FreqResp = frequency_response(hamming_window, window_len)
freq, Voibis_FreqResp = frequency_response(Voibis_windows, window_len)
freq, blackman_FreqResp = frequency_response(blackman_window, window_len)
freq, bartlett_FreqResp = frequency_response(bartlett_window, window_len)
freq, kaiser_FreqRespw = frequency_response(kaiser_window, window_len)
plt.figure()
plt.title("Frequency response")
plt.plot(freq, Rectangle_FreqResp, label="Rectangle")
plt.plot(freq, hanning_FreqResp, label="hanning")
plt.plot(freq, sqrt_hanning_FreqResp, label="sqrt_hanning")
plt.plot(freq, hamming_FreqResp, label="hamming")
plt.plot(freq, Voibis_FreqResp, label="Voibis")
plt.plot(freq, blackman_FreqResp, label="blackman")
plt.plot(freq, bartlett_FreqResp, label="bartlett")
plt.plot(freq, kaiser_FreqRespw, label="kaiser")
plt.ylabel("Magnitude [dB]")
plt.xlabel("Normalized frequency [cycles per sample]")
plt.legend()
plt.tight_layout()
plt.show()
绘制 窗函数和对应的频率响应
1 如何选择窗函数
窗函数频谱的主瓣尽量窄，能量尽可能集中在主瓣内，在频谱分析时能获得较高的频率分辨率
旁瓣增益小且随衰减快，以减小频谱分析时的泄漏失真
但主瓣既窄，旁辨又小衰减又快的窗函数是不容易找到的，比如矩形窗的主瓣宽度最窄，但旁瓣很大，因此在分析处理对应数据时，需要做综合考虑。
下图为针对特定的一段语音信号，加矩形窗与汉宁窗的时域波形及频谱图，Fs=8kHz，窗长取256。可以看出，采用矩形窗时，基音谐波的各个峰都比较尖锐，且整个频谱图显得比较破碎，这是因为矩形窗的主瓣较窄，具有较高的频率分辨率，但是其旁瓣增益较高，因而使基音的相邻皆波之间的干扰比较严重。在相邻谐波间隔内有时叠加，有时抵消，出现了一种随机变化的现象，相邻谐波之间发生频率泄露和混叠，而相对来说，Hamming窗会好多。
2 周期窗和对称窗
在 MATLAB 中，每一个窗函数都可以选择 ‘symmetric’ 或 ‘periodic’ 类型。
’symmetric’ 类型表示窗函数是对称的，主要用于滤波器的设计。
’periodic’ 类型表示窗函数是周期性的，主要用于频谱分析。
下图分别画出了周期窗和对称窗，蓝色的是周期窗(periodic)，红色的是对称窗(symmetric)。在图形上最大的区别是 对称窗有两个最大值，周期窗的最大值在中间。注意如果做stft的时候使用对称的窗函数是不能完美重建的，会有一个比较小的误差。
下图是8个点的频率响应漏，从图中可以看出， periodic拥有稍微窄一点的主瓣，稍微高一点的旁瓣，和稍微低一点的噪声带宽。
窗长的选择
上面已经说过，帧长一般为10~30ms之间，接下来就具体验证帧长会产生什么影响，为了验证该问题，我们人工造一段很简单的数据进行观察，假设overlap为窗长一半，FFT点数与窗长一致，避免引入补零等情况，即为：
通过上图可以验证：长窗具有较高的频率分辨率，较低的时间分辨率。长窗起到了时间上的平均作用。窗宽的选择需折中考虑。短窗具有较好的时间分辨率，能够提取出语音信号中的短时变化(这常常是分析的目的)，损失了频率分辨率。
在python中有很多库都可以创建窗函数，我们一起来探索一下他们是对称窗还是周期窗(非对称)
numpy的hanning函数是对称的
scipy有hanning函数有sym参数设置，默认是对称的
torch的hanning函数有periodic参数设置，默认是非对称的
# -*- coding:utf-8 -*-
# Author:凌逆战 | Never# Date: 2024/3/8
"""
对比不同库中hann窗函数的实现
如果对称(sym=True)的话，有两个最大值，如果不对称(sym=False)的话，有一个最大值
- numpy的hanning函数是对称的 - scipy有hanning函数有sym参数设置，默认是对称的
- torch的hanning函数有periodic参数设置，默认是非对称的
"""
import numpy as np
import torch
import scipy.signal as signal
window_len = 512
def hann_sym(window_len):
"""对称hann窗"""
win = np.zeros(window_len)
for i in range(window_len):
win[i] = 0.5 - 0.5 * np.cos(2 * np.pi * i / (window_len - 1))
return win
def hann_asym(window_len):
"""非对称hann"""
p_win = np.zeros(window_len)
for i in range(window_len):
p_win[i] = np.sin(np.pi * i / window_len)
p_win[i] = p_win[i] * p_win[i]
return p_win
def my_hann_aysm(win_len):
haning_window = np.hanning(win_len + 1) # 对称的hann窗
out = haning_window[0:-1].astype('float32') # 舍弃最后一个元素
return out
scipy_sym = signal.windows.hann(window_len, sym=True) # 对称的hann窗
scipy_Asym = signal.windows.hann(window_len, sym=False) # 非对称的hann窗
hann_sym_c = hann_sym(window_len)
hann_asym_c = hann_asym(window_len)
my_hann= my_hann_aysm(window_len)
print(np.allclose(scipy_sym, hann_sym_c)) # True
print(np.allclose(scipy_Asym, hann_asym_c)) # True
print(np.allclose(my_hann, hann_asym_c)) # True
numpy_window = np.hanning(window_len) # 说明numpy的hanning函数是对称的
print(np.allclose(numpy_window, scipy_sym)) # True
torch_window = torch.hann_window(window_len) # 非对称
torch_window_periodic = torch.hann_window(window_len, periodic=False) # 非周期=对称
# print(torch.argmax(window_torch))
# 判断两个窗函数是否相等
print(np.allclose(scipy_Asym, torch_window.numpy(), rtol=1e-3)) # True
print(np.allclose(scipy_sym, torch_window_periodic.numpy(), rtol=1e-3)) # True
对比不同库中hann窗函数的实现
3 低延迟非对称窗
这里讲的低延迟非对称窗并不是上文的非对称窗(周期窗)，而是真正图形上的非对称窗。
在STFT中，通常会使用重叠的窗来处理信号，以提高频谱分辨率和减少频谱泄漏。重叠的窗会导致相邻窗之间存在重叠部分，这就需要使用OLA技术来将这些重叠部分合并起来，以恢复原始信号。
在进行重叠相加的过程中，会引入一定的延迟，这是因为在重叠部分的处理过程中，需要考虑到前一个窗口和后一个窗口之间的重叠，以确保信号能够完美重建。因此，延迟的产生主要是由于重叠窗口的处理过程中所引入的时间偏移。因此延迟产生的主要因素就有窗长、重叠比例、以及窗的形状。
算法处理延迟一般是由于OLA决定的，比如一个窗长为512，帧移为256的hann窗，一般在做OLA的时候，在256个点之后，第一个完美重建的点才会出来，因此延迟等于帧移。如果我们想要将算法延迟压缩到32个点(2ms)，第一种方法是使用窗长为64，帧移为32个点的窗，这样我们NFFT=64，会导致频率分辨率很低。第二种方法就是使用低延迟非对称窗。在助听器研究中常使用非对称窗函数。
下面举个例子，sqrthann非对称窗，窗长为512
图2：具有高时间（窗口1）和高频谱分辨率（窗口2）的分析和合成窗，用于窗长为K = 512，M = 64和d=64
延迟等于2M-hop_size，如果M=hop_size，如果延迟等于hop_size。
目前非对称窗窗形状有：Orka窗、Tukey 窗、Asqrt hann 窗
def Orka_forward_window(N1=64, N2=448, hop_size=64, NFFT=512):
analysisWindow = np.zeros(NFFT)
for n in range(NFFT):
if n < N1:
analysisWindow[n] = np.sin(n * np.pi / (2 * N1)) ** 2
elif N1 <= n <= N2:
analysisWindow[n] = 1
elif N2 < n <= N2 + hop_size:
analysisWindow[n] = np.sin(np.pi * (N2 + hop_size - n) / (2 * hop_size))
return analysisWindow
def Orka_backward_window(N1=64, N2=448, hop_size=64, NFFT=512):
synthesisWindow = np.zeros(NFFT)
for n in range(NFFT):
if n < N2 - hop_size:
synthesisWindow[n] = 0
elif N2 - hop_size <= n <= N2:
synthesisWindow[n] = np.cos(np.pi * (n - N2) / (2 * hop_size)) ** 2
elif N2 < n <= N2 + hop_size:
synthesisWindow[n] = np.sin(np.pi * (N2 + hop_size - n) / (2 * hop_size))
return synthesisWindow
Orka窗
def TukeyAW(n, N, alpha):
# assert n >= 0
if n < alpha * N:
return 0.5 * (1 - np.cos(np.pi * n / (alpha * N)))
elif n <= N - alpha * N:
return 1
elif n <= N:
return 0.5 * (1 - np.cos(np.pi * (N - n) / (alpha * N)))
def getTukeyAnalysisWindow(filter_length, alpha):
analysisWindow = np.zeros(filter_length)
for i in range(filter_length):
analysisWindow[i] = TukeyAW(i, filter_length, alpha)
return analysisWindow
def getTukeySynthesisWindow(N, A, B, alpha):
synthesisWindow = np.zeros(A)
for i in range(A):
x = N - A + i
numerator = TukeyAW(x, N, alpha)
denonminator = 0
for k in range(int(A / B)):
y = N - A + i % B + k * B
denonminator += TukeyAW(y, N, alpha) ** 2
synthesisWindow[i] = numerator / denonminator
synthesisWindow = np.pad(synthesisWindow, (N - A, 0), 'constant', constant_values=0)
return synthesisWindow
Tukey窗
def getAsqrtAnalysisWindow(N, M, d):
# filter_length, hop_length, d
risingSqrtHann = np.sqrt(np.hanning(2 * (N - M - d) + 1)[:(N - M - d)])
fallingSqrtHann = np.sqrt(np.hanning(2 * M + 1)[:2 * M]) # 下降
window = np.zeros(N)
window[:d] = 0
window[d:N - M] = risingSqrtHann[:N - M - d]
window[N - M:] = fallingSqrtHann[-M:]
return window
def getAsqrtSynthesisWindow(N, M, d):
risingSqrtHannAnalysis = np.sqrt(np.hanning(2 * (N - M - d) + 1)[:(N - M - d)])
fallingSqrtHann = np.sqrt(np.hanning(2 * M + 1)[:2 * M])
risingNoramlizedHann = np.hanning(2 * M + 1)[:M] / risingSqrtHannAnalysis[N - 2 * M - d:N - M - d]
window = np.zeros(N)
window[:-2 * M] = 0
window[-2 * M:-M] = risingNoramlizedHann
window[-M:] = fallingSqrtHann[-M:]
return window
Asqrthann窗
通过OLA过程发现，使用非对称窗确实是在hop_size处信号完美重建，代码见仓库。
参考
【论文】CEC2 E008 Technical Pape【论文】Wang Z Q, Wichern G, Watanabe S, et al. STFT-domain neural speech enhancement with very low algorithmic latency[J]. IEEE/ACM Transactions on Audio, Speech, and Language Processing, 2022, 31: 397-410.【论文】Mauler D, Martin R. A low delay, variable resolution, perfect reconstruction spectral analysis-synthesis system for speech enhancement[C]//2007 15th European Signal Processing Conference. IEEE, 2007: 222-226.