QMF分析综合滤波器

发表于 更新于

分析综合滤波器组的作用

由于音频信号在不同的频率范围具有不同的特性，因此在音频处理之前通常都会使用分析综合滤波器组将音频信号分成不同的频率子带，再分别进行处理。比如，音频编码中常用到的子带编码（Subband coding）；webrtc的VAD中，会对不同的频率子带计算特征，再综合这些特征计算语音概率。

分析综合滤波器组是如何实现

在WebRTC中使用最多的是基于IIR实现的二通道QMF分析综合滤波器组，通过二通道QMF滤波器组，可以很方便实现N等分的子带分解，因此这里只介绍二通道QMF滤波器组。

从上图可以看出，分析综合滤波器包含分析部分和综合部分，当中间没有任何处理时，整个系统的输入输出关系如下：

可以看到分析滤波器部分的高通和低通部分的频率响应网线正好是相对pi/2镜像对称的，QMF叫镜像滤波器的由来。图中Xa0代表的就是输入信号的低频部分，而Xa1代表的就是输入信号的高频部分，这样通过分析滤波器后，就可以对信号的低频部分和高频部分进行分别处理了。

因此只要按下面的等式进行滤波器设计，就可以让A(z)=0，即消除混叠，实现完善重构。

为了效率，通常会采用多相形式实现QMF组，如下图所示，信号处理前都会进行抽取操作，这些实际处理的数据量就减少了，从而提升了执行效率

在QMF组的多相形式中对应的低通滤波器和高通滤波器如上式所示。
WebRTC中的实现

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
void WebRtcSpl_AnalysisQMF(const int16_t* in_data, size_t in_data_length,
                           int16_t* low_band, int16_t* high_band,
                           int32_t* filter_state1, int32_t* filter_state2)
{
    size_t i;
    int16_t k;
    int32_t tmp;
    int32_t half_in1[kMaxBandFrameLength];
    int32_t half_in2[kMaxBandFrameLength];
    int32_t filter1[kMaxBandFrameLength];
    int32_t filter2[kMaxBandFrameLength];
    const size_t band_length = in_data_length / 2;
    RTC_DCHECK_EQ(0, in_data_length % 2);
    RTC_DCHECK_LE(band_length, kMaxBandFrameLength);

    // Split even and odd samples. Also shift them to Q10.
    for (i = 0, k = 0; i < band_length; i++, k += 2)
    {
        half_in2[i] = ((int32_t)in_data[k]) * (1 << 10);
        half_in1[i] = ((int32_t)in_data[k + 1]) * (1 << 10);
    }

    // All pass filter even and odd samples, independently.
    WebRtcSpl_AllPassQMF(half_in1, band_length, filter1,
                         WebRtcSpl_kAllPassFilter1, filter_state1);
    WebRtcSpl_AllPassQMF(half_in2, band_length, filter2,
                         WebRtcSpl_kAllPassFilter2, filter_state2);

    // Take the sum and difference of filtered version of odd and even
    // branches to get upper & lower band.
    for (i = 0; i < band_length; i++)
    {
        tmp = (filter1[i] + filter2[i] + 1024) >> 11;
        low_band[i] = WebRtcSpl_SatW32ToW16(tmp);

        tmp = (filter1[i] - filter2[i] + 1024) >> 11;
        high_band[i] = WebRtcSpl_SatW32ToW16(tmp);
    }
}

上面是WebRTC中关于分析滤波器部分的实现，从代码中可以看出WebRTC中的分析综合滤波器是基于全通滤波器的QMF多相实现，其中的全通滤器采用了IIR实现，即其中的P0(z)和P1(z)都是全通滤波器。通过参考[2]我们可以梳理这个问题的处理流程，通过分析QMF分析综合滤波器满足完美重构的条件，可以得到H0、H1、G0、G1之间的关系，同时H0和H1是基于pi/2，因此只需要知道H0，最终转化成H0低通滤波器的设计问题。进一步的由于采用了基于IIR的全通滤波器，因此只需要考虑相位失真问题，最终QMF分析综合滤波器问题转换成了滤波器的相位均衡问题。虽然我们知道了设计QMF分析综合滤波器的原理和思路，但是想设计一个完全可用的滤波器还是很有难度的，下面我们直接看下WebRTC中QMF分析综合滤波器的效果，如下图所示，可以看到对应的低通滤波器和高通滤波器都有很窄的过渡带，整个系统的幅值响应几乎接近0dB的，同时除了pi/2附近的频带都是近似线性相位的。

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
import scipy.signal as signal
import numpy as np
import matplotlib.pyplot as plt
import control

def analysis_synthesis_filter():
    filter1_coef = [6418.0 / 65536.0, 36982.0 / 65536.0, 57261.0 / 65536.0]
    filter2_coef = [21333.0 / 65536.0, 49062.0 / 65536.0, 63010.0 / 65536.0]
    ha0_0_b = [1, 0, filter1_coef[0]]
    ha0_0_a = [filter1_coef[0], 0, 1]
    ha0_1_b = [1, 0, filter1_coef[1]]
    ha0_1_a = [filter1_coef[1], 0, 1]
    ha0_2_b = [1, 0, filter1_coef[2]]
    ha0_2_a = [filter1_coef[2], 0, 1]

    ha1_0_b = [1, 0, filter2_coef[0]]
    ha1_0_a = [filter2_coef[0], 0, 1]
    ha1_1_b = [1, 0, filter2_coef[1]]
    ha1_1_a = [filter2_coef[1], 0, 1]
    ha1_2_b = [1, 0, filter2_coef[2]]
    ha1_2_a = [filter2_coef[2], 0, 1]

    ha0_b = np.convolve(ha0_2_b, np.convolve(ha0_1_b, ha0_0_b))
    ha0_a = np.convolve(ha0_2_a, np.convolve(ha0_1_a, ha0_0_a))
    ha1_b = np.convolve([1, 0], np.convolve(ha1_2_b, np.convolve(ha1_1_b, ha1_0_b)))
    ha1_a = np.convolve(ha1_2_a, np.convolve(ha1_1_a, ha1_0_a))

    ha0_sys = control.TransferFunction(ha0_b, ha0_a)
    ha1_sys = control.TransferFunction(ha1_b, ha1_a)
    b0_sys = ha1_sys
    b1_sys = ha0_sys
    delay_sys = control.TransferFunction([1, 0], [1])

    lowpass_sys = 0.5 * (ha0_sys + ha1_sys)
    highpass_sys = 0.5 * (ha0_sys - ha1_sys)
    t_sys = delay_sys / 2 * (ha0_sys * b0_sys + ha1_sys * b1_sys)
    # print(t_sys)
    lowpass_num = lowpass_sys.num[0][0]
    lowpass_den = lowpass_sys.den[0][0]
    highpass_num = highpass_sys.num[0][0]
    highpass_den = highpass_sys.den[0][0]
    t_num = t_sys.num[0][0]
    t_den = t_sys.den[0][0]
    lowpass_w, lowpass_h = signal.freqz(lowpass_num, lowpass_den)
    highpass_w, highpass_h = signal.freqz(highpass_num, highpass_den)
    t_w, t_h = signal.freqz(t_num, t_den)
    fig, axes = plt.subplots(3, 2)
    axes[0, 0].plot(lowpass_w/np.pi, 20*np.log10(np.abs(lowpass_h)))
    axes[0, 1].plot(lowpass_w/np.pi, np.unwrap(np.angle(lowpass_h, deg=True)))
    axes[1, 0].plot(highpass_w/np.pi, 20*np.log10(np.abs(highpass_h)))
    axes[1, 1].plot(highpass_w/np.pi, np.unwrap(np.angle(highpass_h, deg=True)))
    axes[2, 0].plot(t_w/np.pi, 20*np.log10(np.abs(t_h)))
    axes[2, 1].plot(t_w/np.pi, np.unwrap(np.angle(t_h, deg=True)))
    plt.show()

最近的一些心得

最后说下近期的两点小心得：

  1. 不懂的知识点请尽早弄懂它。其实对于QMF分析综合滤波器组，在最开始学习音频编码时就遇到了，但是当时没有深入的去搞清楚，最后还是没有躲过去。所以最近也是恶补了很多基础知识，才大致了解了QMF的设计思路。
  2. 不断地输出也许是应对焦虑的一种方法。人到中年难免焦虑，就不停地去学习去吸收，但是往往又是很低效的，时间花了，却什么也没有留下。特别是现在这个信息爆炸的时代，各种信息流，碎片化阅读，让我们看起来收获了很多，其实什么也没有。其实我们要对自己是几平米的房子有基本的认知，小房子就应该放少的、小的东西，定期对房间进行整理，永远保持一定的空间的留白，才会有喘息和美的余地。

参考

  1. Book:《数字信号处理：理论、算法与应用》
  2. Paper: 具有良好重建特性的正交镜像IIR滤波器组的设计新方法
  3. Paper: IIR QMF-bank design for speech and audio subband coding
  4. Blog: WebRTC VAD 中所用滤波器之分析_book_bbyuan的博客-CSDN博客