C++ Butterworth N阶滤波器设计

Question

我正在寻找一个计算Butterworth Nth滤波器设计系数的函数,如Matlab函数:

[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,窗口sinc,分析滤波器......

该库是跨平台的,即在Linux,MacOS和Windows下运行.从 官方网站:

rtfilter应该在任何POSIX平台(GNU/Linux,MacOSX,BSD ...)和Windows平台上编译和运行.

你可以像这样安装它:

$ sudo aptitude install librtfilter-dev librtfilter1

-dev安装软件包后,您甚至可以找到一个示例(使用Butterworth过滤器用法)/usr/share/doc/librtfilter1/examples/butterworth.c.这个例子(以及相应的Makefile)也可以在这里找到.

特别是你对rtf_create_butterworth()功能感兴趣.您可以通过以下命令访问此功能的文档:

$ man rtf_create_butterworth

或者你可以阅读它在这里.

您可以指定任何过滤器顺序将其作为num_pole参数传递给rtf_create_butterworth()函数(据我记得极数,它与过滤器顺序相同).

UPDATE

该库不提供用于系数计算的外部API.它仅提供实际的过滤功能,因此您可以使用rtf_filter()过滤后获取样本(数据).

但是,您可以在库源中找到系数计算的代码.请参阅compute_cheby_iir()函数.这个函数是static,所以它只能在库本身内部使用.但是,您可以将此功能代码原样复制到项目中并使用它.另外,不要让这个函数的名称混淆你:它与Butterworth滤波器和Chebyshev滤波器系数计算的算法相同.

比方说,你准备了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指南,第33章.