C++ Butterworth N阶滤波器设计

Question

我正在寻找一个像 Matlab 函数一样计算 Butterworth Nth 滤波器设计系数的函数：

[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 is a library that provides a set of routines implementing realtime digital filter for multichannel signals (i.e. filtering multiple signals with the same filter parameters). It implements FIR, IIR filters and downsampler for float and double data type (both for real and complex valued signal). Additional functions are also provided to design few usual filters: Butterworth, Chebyshev, windowed sinc, analytical filter...

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

rtfilter should compile and run on any POSIX platform (GNU/Linux, MacOSX, BSD...) and Windows platforms.

你可以这样安装：

$ 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。

C++ Butterworth N阶滤波器设计

C++ Butterworth Nth order filter design

c++

filtering

signal-processing

更新