c++ - 巴特沃斯 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(分母)的 vector 。

还要求该方法在不同平台下工作 - 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 .