• 语音信号处理中的“窗函数”


    文章代码仓库:https://gitee.com/LXP-Never/window_fun

    窗函数贯穿整个语音信号处理,语音信号是一个非平稳的时变信号,但“**短时间内可以认为语音信号是平稳时不变的,一般 10~30ms**。

    对连续的语音分帧做STFT处理,等价于截取一段时间信号,对其进行周期性延拓,从而变成无限长序列,并对该无限长序列做FFT变换,这一截断并不符合傅里叶变换的定义。因此,会导致频谱泄漏和混叠

    • 频谱泄漏:如果不加窗,默认就是矩形窗,时域的乘积就是频域的卷积,使得频谱以实际频率值为中心, 以窗函数频谱波形的形状向两侧扩散,指某一频点能量扩散到相邻频点的现象,会导致幅度较小的频点淹没在幅度较大的频点泄漏分量中
    • 频谱混叠:会在分段拼接处引入虚假的峰值,进而不能获得准确的频谱情况

    加窗的目的让一帧信号的幅度在两端渐变到 0,渐变对傅里叶变换有好处,可以让频谱上的各个峰更细,不容易糊在一起,从而减轻频谱泄漏和混叠的影响

    加窗的代价一帧信号两端的部分被削弱了,没有像中央的部分那样得到重视。弥补的办法就是相互重叠。相邻两帧的起始位置的时间差叫做帧移,常见的取法是取为帧长的一半

    对于语音,窗函数常选汉宁窗(Hanning)、汉明窗(Hamming)、sqrthann及其改进窗,他们的时域波形和幅频响应如下所示:

    1、汉宁窗(Hann)

    w(n)=0.50.5cos(2πnM1)0nM1

    2、汉明窗(Hamming)

    w(n)=0.540.46cos(2πnM1)0nM1

    # -*- 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 如何选择窗函数

    1. 窗函数频谱的主瓣尽量窄,能量尽可能集中在主瓣内,在频谱分析时能获得较高的频率分辨率
    2. 旁瓣增益小且随衰减快,以减小频谱分析时的泄漏失真

     但主瓣既窄,旁辨又小衰减又快的窗函数是不容易找到的,比如矩形窗的主瓣宽度最窄,但旁瓣很大,因此在分析处理对应数据时,需要做综合考虑。

     下图为针对特定的一段语音信号,加矩形窗与汉宁窗的时域波形及频谱图,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.

     

    作者:凌逆战
    欢迎任何形式的转载,但请务必注明出处。
    限于本人水平,如果文章和代码有表述不当之处,还请不吝赐教。
    本文章不做任何商业用途,仅作为自学所用,文章后面会有参考链接,我可能会复制原作者的话,如果介意,我会修改或者删除。

  • 相关阅读:
    Docker-Dockerfile
    SQLite 中的日期和时间
    深度学习文本预处理利器:Tokenizer详解
    总结ES11—ES13新特性——提升开发效率
    基于强化学习的电动汽车的储能系统优化控制和存储容量优化(Matlab代码实现)
    大数据Java基础——集合(泛型,增强for遍历,静态导入,集合嵌套)
    docker&kubernets中级篇(十)
    ROS2自学笔记:动作
    多商户自营连锁小程序商城的作用是什么
    Vue模板语法(01)
  • 原文地址:https://www.cnblogs.com/LXP-Never/p/18175066