如何使过滤器的结果始终达到最终值？

Question

我有一个用 c 编写的指数滤波器，用于在给定时基和时间常数的情况下将输出值过滤为输入值。

void Exp_Filt(float In, float *Out, float time_base, float time_constant) {
  *Out = In + ((*Out - In) * expf(-time_base / time_constant));
}

我的问题是，如果In==0、Out!=0和time_base明显小于time_constant（两者都>0），那么Out无论调用多少次Exp_Filt都不会达到0.0的值。

例如，下面的循环将永远运行：

for (float Out = 1.0f; Out != 0.0f;)
{
  Exp_Filt(0.0f, &Out, 0.1f, 1.0f);
}

是否有任何方法可以解决 IEEE 754 32 位浮点数的精度舍入误差，该误差使 Out 在不需要大量条件语句的情况下无法达到 0？

Answer 1

因为你正在失去精度，你应该以更高的精度进行中间计算（例如 double）并且只在最后转换回来。

所以，你最好这样做：

void
Exp_Filt(double In, float *Out, double time_base, double time_constant)
{
    *Out = In + ((*Out - In) * exp(-time_base / time_constant));
}

此外，您可能永远无法“完美”收敛到 0.0。所以你的循环可能会更好：

#define ERRVAL 1e-5 // some smallish value ...
for (float Out = 1.0f; Out > ERRVAL;)

更新：

通过 fesetround 设置舍入模式似乎有效果 -- 见下文。特别是，fesetround(FE_DOWNWARD)[一次]确实收敛到零。

此外，根据给定的 -std= 选项，我得到的结果略有不同。

我创建了一个测试程序：

/* all.c */

#include <stdio.h>
#include <math.h>
#include <fenv.h>

#ifndef ITER
#define ITER    1000000
#endif

void
Exp_Filt_orig(float In, float *Out, float time_base, float time_constant)
{
    *Out = In + ((*Out - In) * expf(-time_base / time_constant));
}

void
Exp_Filt_fix1(double In, float *Out, double time_base, double time_constant)
{
    *Out = In + ((*Out - In) * exp(-time_base / time_constant));
}

double
Exp_Filt_fix2(double In, double Out, double time_base, double time_constant)
{
    return In + ((Out - In) * exp(-time_base / time_constant));
}

#define DOLOOP(_fnc,_lim) \
    do { \
        float Out; \
        iter = 0; \
        for (Out = 1.0f;  (_lim) && (iter < ITER); ++iter) \
            _fnc(0.0f, &Out, 0.1f, 1.0f); \
        printf("fnc=%s lim='%s' iter=%d Out=%.9g\n", \
            #_fnc,#_lim,iter,(double) Out); \
    } while (0)

#define DOLOOP2(_fnc,_lim) \
    do { \
        double Out; \
        iter = 0; \
        for (Out = 1.0f;  (_lim) && (iter < ITER); ++iter) \
            Out = _fnc(0.0f, Out, 0.1f, 1.0f); \
        printf("fnc=%s lim='%s' iter=%d Out=%.9g\n", \
            #_fnc,#_lim,iter,(double) Out); \
    } while (0)

#define LIM     1e-5
#define LIM2    1e-10

int
main(void)
{
    int iter;

#ifdef __STDC__VERSION__
    printf("std=%ld\n",__STDC_VERSION__);
#endif
    printf("LIM=%.6g\n",LIM);
    printf("LIM2=%.6g\n",LIM2);
    printf("ITER=%d\n",ITER);
    printf("fegetround=%8.8X\n",fegetround());
#ifdef FE
    fesetround(FE);
    printf("fegetround=%8.8X\n",fegetround());
#endif

    DOLOOP(Exp_Filt_orig,Out != 0.0f);
    DOLOOP(Exp_Filt_orig,Out > LIM);
    DOLOOP(Exp_Filt_fix1,Out > LIM);
    DOLOOP2(Exp_Filt_fix2,Out > LIM);
    DOLOOP2(Exp_Filt_fix2,Out > LIM2);
    DOLOOP2(Exp_Filt_fix2,Out != 0.0);

    return 0;
}

这是一个Makefile：

ALL += allxx all90 all99

ifdef O
  OLVL := -O$(O)
endif

ifdef ITER
  XFLAGS += -DITER=$(ITER)
endif

ifdef FE
  XFLAGS += -DFE=$(FE)
endif

.PHONY: all
all: $(ALL)

allxx: all.c
    cc -o allxx all.c -lm $(OLVL) $(XFLAGS)

all90: all.c
    cc -o all90 all.c -lm $(OLVL) $(XFLAGS) -std=c90

all99: all.c
    cc -o all99 all.c -lm $(OLVL) $(XFLAGS) -std=c99

run: all
    ./all90
    ./all99
    ./allxx

clean:
    rm -f $(ALL)

这是 make clean run 的输出：

rm -f allxx all90 all99
cc -o allxx all.c -lm
cc -o all90 all.c -lm   -std=c90
cc -o all99 all.c -lm   -std=c99
./all90
LIM=1e-05
LIM2=1e-10
ITER=1000000
fegetround=00000000
fnc=Exp_Filt_orig lim='Out != 0.0f' iter=1 Out=0
fnc=Exp_Filt_orig lim='Out > LIM' iter=1 Out=0
fnc=Exp_Filt_fix1 lim='Out > LIM' iter=116 Out=9.16608587e-06
fnc=Exp_Filt_fix2 lim='Out > LIM' iter=116 Out=9.16608615e-06
fnc=Exp_Filt_fix2 lim='Out > LIM2' iter=231 Out=9.28532947e-11
fnc=Exp_Filt_fix2 lim='Out != 0.0' iter=1000000 Out=2.47032823e-323
./all99
LIM=1e-05
LIM2=1e-10
ITER=1000000
fegetround=00000000
fnc=Exp_Filt_orig lim='Out != 0.0f' iter=1000000 Out=7.00649232e-45
fnc=Exp_Filt_orig lim='Out > LIM' iter=116 Out=9.16610225e-06
fnc=Exp_Filt_fix1 lim='Out > LIM' iter=116 Out=9.16608587e-06
fnc=Exp_Filt_fix2 lim='Out > LIM' iter=116 Out=9.16608615e-06
fnc=Exp_Filt_fix2 lim='Out > LIM2' iter=231 Out=9.28532947e-11
fnc=Exp_Filt_fix2 lim='Out != 0.0' iter=1000000 Out=2.47032823e-323
./allxx
LIM=1e-05
LIM2=1e-10
ITER=1000000
fegetround=00000000
fnc=Exp_Filt_orig lim='Out != 0.0f' iter=1000000 Out=7.00649232e-45
fnc=Exp_Filt_orig lim='Out > LIM' iter=116 Out=9.16610225e-06
fnc=Exp_Filt_fix1 lim='Out > LIM' iter=116 Out=9.16608587e-06
fnc=Exp_Filt_fix2 lim='Out > LIM' iter=116 Out=9.16608615e-06
fnc=Exp_Filt_fix2 lim='Out > LIM2' iter=231 Out=9.28532947e-11
fnc=Exp_Filt_fix2 lim='Out != 0.0' iter=1000000 Out=2.47032823e-323

但是，这里是 make FE=FE_DOWNWARD clean run 的输出：

rm -f allxx all90 all99
cc -o allxx all.c -lm  -DFE=FE_DOWNWARD
cc -o all90 all.c -lm  -DFE=FE_DOWNWARD -std=c90
cc -o all99 all.c -lm  -DFE=FE_DOWNWARD -std=c99
./all90
LIM=1e-05
LIM2=1e-10
ITER=1000000
fegetround=00000000
fegetround=00000400
fnc=Exp_Filt_orig lim='Out != 0.0f' iter=1000000 Out=3.40282346e+38
fnc=Exp_Filt_orig lim='Out > LIM' iter=1000000 Out=3.40282346e+38
fnc=Exp_Filt_fix1 lim='Out > LIM' iter=116 Out=9.16604039e-06
fnc=Exp_Filt_fix2 lim='Out > LIM' iter=116 Out=9.16608615e-06
fnc=Exp_Filt_fix2 lim='Out > LIM2' iter=231 Out=9.28532947e-11
fnc=Exp_Filt_fix2 lim='Out != 0.0' iter=7427 Out=0
./all99
LIM=1e-05
LIM2=1e-10
ITER=1000000
fegetround=00000000
fegetround=00000400
fnc=Exp_Filt_orig lim='Out != 0.0f' iter=1015 Out=0
fnc=Exp_Filt_orig lim='Out > LIM' iter=116 Out=9.16599037e-06
fnc=Exp_Filt_fix1 lim='Out > LIM' iter=116 Out=9.16604039e-06
fnc=Exp_Filt_fix2 lim='Out > LIM' iter=116 Out=9.16608615e-06
fnc=Exp_Filt_fix2 lim='Out > LIM2' iter=231 Out=9.28532947e-11
fnc=Exp_Filt_fix2 lim='Out != 0.0' iter=7427 Out=0
./allxx
LIM=1e-05
LIM2=1e-10
ITER=1000000
fegetround=00000000
fegetround=00000400
fnc=Exp_Filt_orig lim='Out != 0.0f' iter=1015 Out=0
fnc=Exp_Filt_orig lim='Out > LIM' iter=116 Out=9.16599037e-06
fnc=Exp_Filt_fix1 lim='Out > LIM' iter=116 Out=9.16604039e-06
fnc=Exp_Filt_fix2 lim='Out > LIM' iter=116 Out=9.16608615e-06
fnc=Exp_Filt_fix2 lim='Out > LIM2' iter=231 Out=9.28532947e-11
fnc=Exp_Filt_fix2 lim='Out != 0.0' iter=7427 Out=0

编辑： 原创 post 正在制作 DV：

当您调用给定函数时，如果函数的 参数之一 是 float，调用函数会自动将其作为 double 传递并被 调用的 函数解释为 double。

也就是说，In、time_base 和 time_constant 是 double，即使您指定了 float.只有 Out [因为它是指向 float 的指针并且不是按值传递的]，是 float.

而且，在表达式中，float 在求值过程中会自动升级为 double。

Answer 2

Is there any way to account for the precision round-off error of IEEE 754 32bit floating-point that is keeping Out from reaching 0 without needing a multitude of conditional statements?

expf(-0.1f / 1.0f) 是 0.904837...，所以 FLT_TRUE_MIN 的任何乘法都会舍入到最近的舍入到 FLT_TRUE_MIN，所以我们陷入了循环。

另一种方法是使用 Monte Carlo rounding

void Exp_Filt_MC(float In, float *Out, float time_base, float time_constant) {
  if (rand()&1) {
    fesetround(FE_UPWARD);
  } else {
    fesetround(FE_DOWNWARD);
  }
  *Out = In + ((*Out - In) * expf(-time_base / time_constant));

  // Note 
  // Better could would restore the rounding mode of function entry.
}

这最终会退出循环，因为有时 FLT_TRUE_MIN*0.904837 会是 0.0.

这里可能有主题的变化，或者限制为仅在 (*Out - In) 低于正常水平时才这样做。

然而“不需要大量的条件语句？”这似乎是最简单的。

我也喜欢这样的想法，如果 *Out 在函数的末尾是次正规的，就是 return 0.0f.

这类似于

for (float Out = 1.0f; fabsf(Out) <= FLT_MIN; )

如何使过滤器的结果始终达到最终值？

How to get result of filter to always go to final value?

c

floating-point