C++ Butterworth N 阶滤波器设计

Question

我正在寻找一个像 Matlab 函数一样计算巴特沃斯第 N 个滤波器设计系数的函数：

[bl,al]=butter(but_order,Ws); 

和

[bh,ah]=butter(but_order,2*bandwidth(1)/fs,'high');

我发现了许多计算二阶而不是 N 阶的示例（例如，我使用 18 阶...）。 - 不幸的是我对DSP一无所知。

您知道任何库或轻松实现此方法的方法吗？当我知道只是命令时，切断频率和采样率。我只需要得到 B（分子）和 A（分母）的向量。

还要求该方法在不同的平台下工作 - Windows、Linux、...

提前致谢。

Answer 1

很容易找到（在 Debian 或 Ubuntu 中）：

$ aptitude search ~dbutterworth | grep lib

马上给你答案：

p   librtfilter-dev         - realtime digital filtering library (dev)
p   librtfilter1            - realtime digital filtering library        
p   librtfilter1-dbg        - realtime digital filtering library (dbg)

所以你需要名为rtfilter 的库。说明：

rtfilter 是一个库，它提供了一组实现多通道信号实时数字滤波器的例程（即用相同的滤波器参数过滤多个信号）。它实现了浮点数和双精度数据类型的 FIR、IIR 滤波器和下采样器（适用于实值和复值信号）。还提供了额外的函数来设计一些常用的滤波器：Butterworth、Chebyshev、windowed sinc、分析滤波器...

这个库是跨平台的，即在 Linux、MacOS 和 Windows 下工作。从 official site:

rtfilter 应该在任何 POSIX 平台（GNU/Linux、MacOSX、BSD...）和 Windows 平台上编译和运行。

你可以这样安装：

$ sudo aptitude install librtfilter-dev librtfilter1

安装-dev 包后，您甚至可以在/usr/share/doc/librtfilter1/examples/butterworth.c 找到一个示例（使用Butterworth 过滤器）。这个例子（连同对应的Makefile）也可以在here找到。

您对rtf_create_butterworth() 函数特别感兴趣。您可以通过命令访问此函数的文档：

$ man rtf_create_butterworth

或者你可以阅读它here。

您可以指定任何过滤顺序将其作为num_pole 参数传递给rtf_create_butterworth() 函数（据我所知，极数与过滤顺序相同）。

更新

此库不提供用于系数计算的外部 API。它只提供实际的过滤能力，所以您可以使用rtf_filter()获取过滤后的样本（数据）。

但是，你可以在库中找到系数计算的代码。请参阅compute_cheby_iir() 函数。这个函数是static，所以它只能在库本身内部使用。但是，您可以按原样将该功能代码复制到您的项目中并使用它。另外，不要让这个函数的名称混淆你：它是巴特沃斯滤波器和切比雪夫滤波器系数计算的相同算法。

假设你已经为rtf_create_butterworth()函数准备了参数：

const double cutoff     = 8.0;          /* cutoff frequency, in Hz */
const double fs         = 512.0;        /* sampling rate, in Hz */

unsigned int nchann     = 1;            /* channels number */
int proctype            = RTF_FLOAT;    /* samples have float type */
double fc               = cutoff / fs;  /* normalized cut-off frequency, Hz */
unsigned int num_pole   = 2;            /* filter order */
int highpass            = 0;            /* lowpass filter */

现在您要计算过滤器的分子和分母。我已经为你编写了包装器：

struct coeff {
    double *num;
    double *den;
};

/* TODO: Insert compute_cheby_iir() function here, from library:
 * https://github.com/nbourdau/rtfilter/blob/master/src/common-filters.c#L250
 */

/* Calculate coefficients for Butterworth filter.
 * coeff: contains calculated coefficients
 * Returns 0 on success or negative value on failure.
 */
static int calc_coeff(unsigned int nchann, int proctype, double fc,
                      unsigned int num_pole, int highpass,
                      struct coeff *coeff)
{
    double *num = NULL, *den = NULL;
    double ripple = 0.0;
    int res = 0;

    if (num_pole % 2 != 0)
        return -1;

    num = calloc(num_pole+1, sizeof(*num));
    if (num == NULL)
        return -2;
    den = calloc(num_pole+1, sizeof(*den));
    if (den == NULL) {
        res = -3;
        goto err1;
    }

    /* Prepare the z-transform of the filter */
    if (!compute_cheby_iir(num, den, num_pole, highpass, ripple, fc)) {
        res = -4;
        goto err2;
    }

    coeff->num = num;
    coeff->den = den;

    return 0;

err2:
    free(den);
err1:
    free(num);
    return res;
}

你可以像这样使用这个包装器：

int main(void)
{
    struct coeff coeff;
    int res;
    int i;

    /* Calculate coefficients */
    res = calc_coeff(nchann, proctype, fc, num_pole, highpass, &coeff);
    if (res != 0) {
        fprintf(stderr, "Error: unable to calculate coefficients: %d\n", res);
        return EXIT_FAILURE;
    }

    /* TODO: Work with calculated coefficients here (coeff.num, coeff.den) */
    for (i = 0; i < num_pole+1; ++i)
        printf("num[%d] = %f\n", i, coeff.num[i]);
    for (i = 0; i < num_pole+1; ++i)
        printf("den[%d] = %f\n", i, coeff.den[i]);

    /* Don't forget to free memory allocated in calc_coeff() */
    free(coeff.num);
    free(coeff.den);

    return EXIT_SUCCESS;
}

如果您对这些系数计算的数学背景感兴趣，请查看DSP Guide, chapter 33。