28、高性能数字滤波器：整数滤波器的设计与应用

高性能数字滤波器：整数滤波器的设计与应用

1. 引言

在数字信号处理中，滤波器是一种重要的工具，用于去除信号中的噪声或提取特定频率的信号。传统的浮点滤波器在一些微控制器平台上可能效率较低，因为许多小型微控制器缺乏硬件浮点单元（FPU），浮点运算需要软件模拟，这会导致效率低下。因此，为了在嵌入式系统中实现高性能滤波，我们需要设计使用整数运算的滤波器。

2. 一阶整数滤波器

2.1 滤波器表达式

在实现整数滤波器时，首先需要将浮点样本值和系数表示为归一化的整数值。一阶 FIR 滤波器的表达式可以重写为整数形式：
[y_1 = \frac{\beta_0 x_0 + \beta_1 x_1 + \frac{1}{2}(\beta_0 + \beta_1)}{\beta_0 + \beta_1}]
其中，(y_1)、(x_0)、(x_1)、(\beta_0) 和 (\beta_1) 是无符号整数值，分子中的额外项 (\frac{1}{2}(\beta_0 + \beta_1)) 用于处理无符号整数的舍入。

2.2 模板类实现

以下是一个可扩展的模板类，用于实现上述一阶整数滤波器：

template<const std::size_t resol = 4U,
         typename sample_t = std::uint16_t,
         typename value_t = sample_t,
         typename result_t = sample_t>
class fir_01
{
public:
    typedef sample_t sample_type;
    typedef value_t value_type;
    typedef result_t result_type;
    typedef std::int_fast16_t weight_type;

    fir_01(const sample_type& val = 0U)
        : result(val * resol)
    {
        std::fill(values.begin(), values.end(), result);
    }

    template<const weight_type B0, const weight_type B1>
    void new_sample(const sample_type& val)
    {
        values[0U] = values[1U];
        values[1U] = val * static_cast<value_type>(resol);
        value_type new_val = (B0 * values[0U]) + (B1 * values[1U]);
        result = (new_val + ((B0 + B1) / 2)) / (B0 + B1);
    }

    result_type get_result() const
    {
        return (result + (resol / 2U)) / resol;
    }

private:
    result_type result;
    std::array<value_type, 2U> values;
};

2.3 模板参数说明

• resol ：分辨率缩放参数，通过将每个新样本乘以一个常量整数来提供分辨率缩放。较高的 resol 值可以更接近模拟滤波器的效果。在 get_result() 函数中，会从滤波器结果中移除分辨率缩放。
• sample_t 、 value_t 和 result_t ：缩放参数，用于定义滤波器的维度。它们可以设置为 8 位、16 位、32 位甚至 64 位，提供了多个自由度的可扩展性，因为表示滤波器样本、延迟线和滤波器结果的变量大小可以独立设置。

2.4 使用示例

typedef fir_01<> filter_type;
filter_type f(4U);

void do_something()
{
    // The result of the first call is 16.
    f.new_sample<7, 1>(100U);
}

在这个示例中，创建了一个 fir_01 对象 f ，并使用初始值 4 进行初始化。调用 new_sample<7, 1>(100U) 时，新样本值 100 被放置在延迟线的顶部，并乘以系数 1，旧值 4 乘以系数 7。最终结果为：
[\frac{(7 \times 4) + (1 \times 100) + (8 / 2)}{8} = 16]

2.5 性能优化

编译器可以对 fir_01 类的代码进行高度优化。例如，构造函数中的 std::fill() 循环可以被展开，乘法运算可以被替换为快速的移位和加法操作，除法运算可以被替换为右移操作。这使得滤波器代码非常高效，接近汇编编程的效率。

3. N 阶整数 FIR 滤波器

3.1 滤波器定义

N 阶 FIR 滤波器由差分方程定义：
[y_n = b_0 x[n] + b_1 x[n - 1] + \cdots + b_N x[n - N] = \sum_{i = 0}^{N} b_i x[n - i]]
其中，(x[n]) 是延迟线的值，(y_n) 是滤波器结果，(b_i) 是系数，(N) 是滤波器阶数。

3.2 模板类实现

以下是一个用于实现 N 阶整数 FIR 滤波器的模板类：

template<const std::size_t order,
         const std::size_t resol = 4U,
         typename sample_t = std::uint16_t,
         typename value_t = sample_t,
         typename result_t = sample_t>
class fir_order_n
{
public:
    static_assert((order > 0U) && (order < 48U),
                  "error: filter order must be from 1 to 48");

    fir_order_n() { }
    explicit fir_order_n(const sample_t&) { }

    template<typename... dummy_parameters>
    void new_sample(const sample_t&) { }

    result_t get_result() const { return result_t(0); }
};

3.3 模板类特化

由于每个不同的滤波器阶数 N 需要有自己特定的 new_sample() 子例程模板变体，因此每个 N 阶 fir_order_n 类都需要显式特化。以下是一个 5 阶滤波器的特化示例：

template<const std::size_t resol,
         typename sample_t,
         typename value_t,
         typename result_t>
class fir_order_n<5U, resol, sample_t, value_t, result_t>
{
public:
    typedef sample_t sample_type;
    typedef value_t value_type;
    typedef result_t result_type;

    fir_order_n() : result(0)
    {
        std::fill(data.begin(), data.end(), result);
    }

    explicit fir_order_n(const sample_type& x)
        : result(value_type(x) * resol)
    {
        std::fill(data.begin(), data.end(), result);
    }

    template<const std::int_fast16_t B0,
             const std::int_fast16_t B1,
             const std::int_fast16_t B2,
             const std::int_fast16_t B3,
             const std::int_fast16_t B4,
             const std::int_fast16_t B5>
    void new_sample(const sample_type& x)
    {
        // Shift the delay line.
        std::copy(data.begin() + 1U, data.end(), data.begin());
        // Store the new sample at top of delay line.
        *(data.end() - 1U) = value_type(x) * resol;
        // Calculate the FIR algorithm.
        const value_type new_val =
            value_type(data[0U] * B0) +
            value_type(data[1U] * B1) +
            value_type(data[2U] * B2) +
            value_type(data[3U] * B3) +
            value_type(data[4U] * B4) +
            value_type(data[5U] * B5);
        constexpr std::int_fast16_t weight = B0 + B1 + B2 + B3 + B4 + B5;
        result = (new_val + (weight / 2)) / weight;
    }

    result_type get_result() const
    {
        return (result + (resol / 2U)) / resol;
    }

private:
    result_type result;
    std::array<value_type, 6U> data;
};

3.4 使用示例

typedef fir_order_n<5U> filter_type;
filter_type f(4U);

void do_something()
{
    f.new_sample<5, 5, 6, 6, 5, 5>(100U);
}

在这个示例中，使用一个 5 阶低通滤波器对样本值 100 进行滤波。滤波器的系数为 5、5、6、6、5、5，系数和为 32，分辨率缩放为 4。最终结果为：
[\frac{(5 \times 4) + (5 \times 4) + (6 \times 4) + (6 \times 4) + (5 \times 4) + (5 \times 100) + 16}{32} = 19]

3.5 性能分析

与一阶滤波器类似，N 阶滤波器的 new_sample() 函数也可以被编译器高度优化。然而，高阶滤波器通常需要更多的代码和运行时间，特别是当滤波器的维度较大时。例如，一个 17 阶的滤波器可能需要大量的 32 位操作，这对于 8 位或 16 位的目标平台来说可能过于复杂。在这种情况下，可以考虑使用多个低阶滤波器级联来实现相同的滤波效果，从而减少 CPU 负载。

4. 滤波器设计流程

以下是设计一个低通 FIR 滤波器的流程图：

graph TD;
    A[确定信号特性] --> B[选择采样频率];
    B --> C[确定滤波器参数];
    C --> D[使用滤波器设计工具计算系数];
    D --> E[选择滤波器阶数和维度];
    E --> F[实现滤波器代码];
    F --> G[测试滤波器性能];
    G --> H{性能是否满足要求};
    H -- 是 --> I[完成设计];
    H -- 否 --> C;

4.1 确定信号特性

首先需要了解待滤波信号的主要频率成分、噪声频率和幅度等特性。例如，一个信号的主要成分是 60 Hz 的正弦波，叠加了 600 Hz 的噪声。

4.2 选择采样频率

根据信号的噪声频率，选择合适的采样频率。通常，建议每半个噪声周期采样 3 - 4 个点。例如，对于 600 Hz 的噪声，采样频率可以选择为 4000 Hz。

4.3 确定滤波器参数

选择滤波器的通带和阻带频率、阻带衰减和通带波纹等参数。例如，选择通带上限频率为 200 Hz，阻带下限频率为 600 Hz，阻带衰减为 40 dB，通带波纹为 1 dB。

4.4 使用滤波器设计工具计算系数

将滤波器参数输入到滤波器设计工具中，计算出滤波器的系数。例如，对于一个 18 阶的滤波器，工具会计算出 18 个双精度系数。

4.5 选择滤波器阶数和维度

根据滤波器的性能要求和目标平台的硬件特性，选择合适的滤波器阶数和维度。例如，对于 16 位的微控制器，可以选择 16 位的无符号整数类型来表示样本和结果。

4.6 实现滤波器代码

根据选择的滤波器阶数和维度，实现滤波器代码。可以使用前面介绍的模板类来实现整数滤波器。

4.7 测试滤波器性能

使用测试数据对滤波器进行测试，评估滤波器的性能。例如，使用一个包含 101 个样本的测试数据集，对滤波器进行测试。

4.8 调整设计

如果滤波器的性能不满足要求，需要返回前面的步骤，调整滤波器参数或阶数，直到性能满足要求为止。

5. 总结

本文介绍了一阶和 N 阶整数 FIR 滤波器的设计和实现方法。通过使用整数运算，这些滤波器可以在嵌入式系统中实现高性能滤波。同时，我们还介绍了滤波器设计的流程，包括确定信号特性、选择采样频率、确定滤波器参数、计算系数、选择滤波器阶数和维度、实现滤波器代码和测试滤波器性能等步骤。在实际应用中，需要根据具体的需求和目标平台的硬件特性，选择合适的滤波器设计和实现方法。

6. 实际滤波器示例

6.1 信号特性

考虑一个未经过滤波的原始信号 S，它可能是由电压测量输入到 10 位 ADC 得到的。该信号的主要成分是一个频率为 60 Hz、振幅为 150、偏移为 250 的正弦波，同时叠加了一个强异步噪声。噪声的频率是信号频率的 10 倍（600 Hz），振幅是信号振幅的 1/5（30），偏移为 0.317 ms。其数学表达式为：
[S = 250 + 150 \times \left(\sin(0.12\pi t) + \frac{1}{5} \sin(0.317 + 1.2\pi t)\right)]
其中，t 是时间（单位：ms）。

6.2 测试数据

为了测试滤波器性能，使用采样频率为 4 kHz 生成了 101 个数字化点作为测试数据，存储在一个静态常量 STL 数组中：

#include <cstdint>
#include <array>

const std::array<std::uint16_t, 101U> data =
{
    {
        250U, 288U, 306U, 301U, 287U, 288U, 312U, 351U,
        381U, 386U, 371U, 354U, 357U, 381U, 412U, 428U,
        417U, 390U, 370U, 372U, 392U, 411U, 409U, 383U,
        347U, 326U, 328U, 343U, 350U, 333U, 296U, 258U,
        241U, 246U, 258U, 256U, 231U, 190U, 158U, 150U,
        162U, 176U, 170U, 141U, 106U,  87U,  93U, 116U,
        132U, 125U, 100U,  77U,  75U,  97U, 129U, 147U,
        141U, 123U, 113U, 127U, 162U, 198U, 215U, 209U,
        195U, 197U, 224U, 264U, 297U, 306U, 296U, 285U,
        293U, 325U, 363U, 386U, 383U, 364U, 352U, 363U,
        392U, 420U, 427U, 409U, 381U, 368U, 377U, 400U,
        414U, 403U, 371U, 338U, 324U, 332U, 348U, 348U,
        322U, 282U, 250U, 240U, 250U
    }
};

6.3 17 阶低通滤波器

6.3.1 滤波器设计

选择通带上限频率为 200 Hz，阻带下限频率为 600 Hz，阻带衰减为 40 dB，通带波纹为 1 dB。使用滤波器设计工具计算得到 18 个双精度系数，对应于 17 阶滤波器的缩放整数系数。

6.3.2 代码实现

#include <iostream>
#include <math/filters/fir_order_n.h>

typedef fir_order_n<17U, 64U, std::uint16_t, std::uint32_t> filter_type;

void do_something()
{
    filter_type f(data[0U]);
    std::cout << f.get_result() << "\n";
    std::for_each(
        data.begin() + 1U,
        data.end(),
        [&f](const std::uint16_t& s)
        {
            f.new_sample<-2, -2, -2, -1, 3, 9, 15, 20, 24,
                          24, 20, 15, 9, 3, -1, -2, -2, -2>(s);
            std::cout << f.get_result() << "\n";
        });
}

6.3.3 性能分析

该 17 阶滤波器在 32 位目标平台上运行速度较快，例如在某 32 位微控制器上， new_sample() 函数大约需要 9.6 μs。由于采样率为 4 kHz，采样周期为 250 μs，滤波器操作大约占用总 CPU 功率的 3.8%。然而，该滤波器有许多 32 位操作，对于 8 位或 16 位架构来说过于复杂。在 8 位目标平台上， new_sample() 函数大约需要 56 μs，占用总 CPU 功率的 22%，CPU 负载过高。

6.4 级联 5 阶低通滤波器

6.4.1 滤波器设计

使用与 17 阶滤波器相同的滤波器参数，但将滤波器阶数限制为 5，得到系数为 (5, 5, 6, 6, 5, 5)。

6.4.2 代码实现

#include <iostream>

typedef fir_order_n<5U> filter_type;

void do_something()
{
    filter_type f1(data[0U]);
    filter_type f2(f1.get_result());
    std::cout << f2.get_result() << std::endl;
    std::for_each(
        data.begin() + 1U,
        data.end(),
        [&f1, &f2](const std::uint16_t& s)
        {
            f1.new_sample<5, 5, 6, 6, 5, 5>(s);
            f2.new_sample<5, 5, 6, 6, 5, 5>(f1.get_result());
            std::cout << f2.get_result() << std::endl;
        });
}

6.4.3 性能优势

通过使用两个 5 阶低通滤波器级联，可以在 8 位或 16 位架构上显著降低 CPU 负载，同时实现与 17 阶滤波器相近的滤波效果。

6.5 不同滤波器性能对比

滤波器类型	阶数	目标平台	运行时间（μs）	CPU 负载（4 kHz 采样率）	滤波效果
17 阶低通滤波器	17	32 位	9.6	3.8%	优秀
17 阶低通滤波器	17	8 位	56	22%	优秀
级联 5 阶低通滤波器	5 + 5	8 位/16 位	较低	较低	接近 17 阶

7. 滤波器使用注意事项

7.1 模板参数选择

• resol ：应选择 2 的倍数，最好是 (2^n)（n 为小正整数），以确保舍入校正精确，并让编译器可以用高效的移位操作替代除法。
• 滤波器系数 (B_i)：系数之和 (\sum_{i=0}^{N} |B_i|) 应是 2 的小整数次幂，以便编译器进行优化。
• 数据类型：根据样本值范围和目标平台硬件特性，选择合适的 sample_t 、 value_t 和 result_t ，避免数值溢出。

7.2 滤波器阶数选择

• 高阶滤波器通常具有更好的滤波性能，但需要更多的计算资源和运行时间。
• 对于资源受限的平台，可以考虑使用多个低阶滤波器级联来实现相同的滤波效果。

7.3 性能优化

• 编译器可以对滤波器代码进行高度优化，如展开循环、替换乘法和除法为移位和加法操作等。在编写代码时，应尽量利用这些优化特性。

8. 总结与展望

8.1 总结

本文详细介绍了一阶和 N 阶整数 FIR 滤波器的设计、实现和应用。通过使用整数运算，这些滤波器可以在嵌入式系统中实现高性能滤波。同时，通过实际滤波器示例展示了不同阶数和架构下滤波器的性能差异，并提供了滤波器设计和使用的注意事项。

8.2 展望

在未来的数字信号处理中，随着硬件技术的不断发展，滤波器的设计和实现将更加高效和灵活。例如，新的处理器架构可能支持更高效的整数运算，从而进一步提高滤波器的性能。此外，深度学习等新兴技术也可能为滤波器设计带来新的思路和方法。

总之，整数 FIR 滤波器在嵌入式系统中具有重要的应用价值，通过合理的设计和优化，可以满足不同应用场景的需求。在实际应用中，需要根据具体的需求和目标平台的硬件特性，选择合适的滤波器设计和实现方法。