20 巴特沃斯滤波器的设计

发布于 2021-10-25 15:26:17 浏览：1228 订阅该版

有没有软件包支持 “C语言巴特沃斯滤波器”？或者哪位大佬用过，给个思路。

myLLgf 2021-10-26

这家伙很懒，什么也没写！

给你个二阶低通巴特沃斯的代码，记得采纳

typedef struct lowpassFilter2p
{
    // 历史采样数据
    float _delay_element1;
    float _delay_element2;
    // 权重系数
    float _a1;
    float _a2;
    float _b0;
    float _b1;
    float _b2;
    // 截至频率
    float _cutoff_freq;
    // 采样频率
    float _sample_freq;
} lowpassFilter2p_t;
void lp2p_init(lowpassFilter2p_t *lp, float sample_freq, float cutoff_freq)
{
    lp2p_set_cutoff_frequency(lp, sample_freq, cutoff_freq);
}
float lp2p_reset(lowpassFilter2p_t *lp, float sample)
{
    const float input = isfinite(sample) ? sample : 0.0f;
    if (fabsf(1 + lp->_a1 + lp->_a2) > FLT_EPSILON)
    {
        lp->_delay_element1 = lp->_delay_element2 = input / (1 + lp->_a1 + lp->_a2);
        if (!isfinite(lp->_delay_element1) || !isfinite(lp->_delay_element2))
        {
            lp->_delay_element1 = lp->_delay_element2 = input;
        }
    }
    else
    {
        lp->_delay_element1 = lp->_delay_element2 = input;
    }
    return lp2p_apply(lp, input);
}
void lp2p_disable(lowpassFilter2p_t *lp)
{
    lp->_sample_freq = 0.f;
    lp->_cutoff_freq = 0.f;
    lp->_delay_element1 = 0.0f;
    lp->_delay_element2 = 0.0f;
    lp->_b0 = 1.f;
    lp->_b1 = 0.f;
    lp->_b2 = 0.f;
    lp->_a1 = 0.f;
    lp->_a2 = 0.f;
}
void lp2p_set_cutoff_frequency(lowpassFilter2p_t *lp, float sample_freq,
                               float cutoff_freq)
{
    if ((sample_freq <= 0.f) || (cutoff_freq <= 0.f) ||
        (cutoff_freq >= sample_freq / 2) || !isfinite(sample_freq) ||
        !isfinite(cutoff_freq))
    {
        lp2p_disable(lp);
        return;
    }
    // reset delay elements on filter change
    lp->_delay_element1 = 0;
    lp->_delay_element2 = 0;
    lp->_cutoff_freq = constrain_float(
        cutoff_freq, 1.f,
        sample_freq / 2); // TODO: min based on actual numerical limit
    lp->_sample_freq = sample_freq;
    const float fr = lp->_sample_freq / lp->_cutoff_freq;
    const float ohm = tanf(M_PI_F / fr);
    const float c = 1.f + 2.f * cosf(M_PI_F / 4.f) * ohm + ohm * ohm;
    lp->_b0 = ohm * ohm / c;
    lp->_b1 = 2.f * lp->_b0;
    lp->_b2 = lp->_b0;
    lp->_a1 = 2.f * (ohm * ohm - 1.f) / c;
    lp->_a2 = (1.f - 2.f * cosf(M_PI_F / 4.f) * ohm + ohm * ohm) / c;
    if (!isfinite(lp->_b0) || !isfinite(lp->_b1) || !isfinite(lp->_b2) ||
        !isfinite(lp->_a1) || !isfinite(lp->_a2))
    {
        lp2p_disable(lp);
    }
}
float lp2p_apply(lowpassFilter2p_t *lp, float sample)
{
    float delay_element_0 =
        sample - lp->_delay_element1 * lp->_a1 - lp->_delay_element2 * lp->_a2;
    const float output = delay_element_0 * lp->_b0 +
                         lp->_delay_element1 * lp->_b1 +
                         lp->_delay_element2 * lp->_b2;
    lp->_delay_element2 = lp->_delay_element1;
    lp->_delay_element1 = delay_element_0;
    return output;
}

yc961213911 2021-10-27

这家伙很懒

大哥，这两处报错怎么处理？

lp->_cutoff_freq = constrain_float(cutoff_freq, 1.f,sample_freq / 2); // TODO: min based on actual numerical limit

这个函数怎么实现？
“M_PI_F”这个是圆周率π吗？

myLLgf 2021-10-27

这家伙很懒，什么也没写！

@yc961213911 一个是圆周率，另一个是限幅

static inline float constrain_float(float x, float low, float high)
{
    if (x < low)
        return low;
    else if (x > high)
        return high;
    else
        return x;
}

yc961213911 2021-10-28

这家伙很懒

大佬，函数怎么用呢？
比如我对一个频率在50HZ~80HZ之间的波形进行采样，采集了160个点放在一个数组。
然后把每个数组元素都调用lp2p_apply()这个函数重新计算出一个值，代替原来的值吗？

myLLgf 2021-10-29

这家伙很懒，什么也没写！

@yc961213911 先 lp2p_init 初始化，然后用 lp2p_apply 依次遍历所有采样点，它的返回值是滤波结果

近我者富 2024-03-13

这家伙很懒，什么也没写！

if (fabsf(1 + lp->_a1 + lp->_a2) > FLT_EPSILON)
    {
        lp->_delay_element1 = lp->_delay_element2 = input / (1 + lp->_a1 + lp->_a2);
        if (!isfinite(lp->_delay_element1) || !isfinite(lp->_delay_element2))
        {
            lp->_delay_element1 = lp->_delay_element2 = input;
        }
    }
    else
    {
        lp->_delay_element1 = lp->_delay_element2 = input;
    }