如何让过滤器的结果始终达到最终值？答案

【问题标题】：How to get result of filter to always go to final value?如何让过滤器的结果始终达到最终值？
【发布时间】：2021-08-09 22:06:09
【问题描述】：

我有一个用 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无论如何都不会达到0.0的值多次调用 Exp_Filt。

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

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

有没有什么方法可以解释 IEEE 754 32 位浮点的精度舍入误差，使 Out 在不需要大量条件语句的情况下达到 0？

【问题讨论】：

How dangerous is it to compare floating point values?
过滤器一般用于建模或控制物理现象。在某些时候，Out 和它的目标In 之间的差异在物理上变得微不足道——声音差异低于人类的分辨能力，无线电信号在噪音中丢失，机器无法如此精确地控制。一般Out != 0.0f可以换成fabsf(Out-In) < threshold。但这出现的背景是什么？如果目标只是达到for (float Out = 1.0f; Out != 0.0f;){ Exp_Filt(0.0f, &Out, 0.1f, 1.0f); } 的最终状态，那么使用Out = In 很容易完成。
经过一些迭代后，您将 pow(2,min_exp) 乘以 1 到 0.5 之间的数字，然后向上舍入到 pow(2,min_exp)（在大多数舍入行为下）。
@chtz：这很明显；问题是如何处理它。此外，“min_exp”通常被认为是最小正常值的指数，如 C 标准的FLT_MIN_EXP（有效数的模缩放）。所讨论的值是最小正值，2^(1-p+最小指数)，其中 p 是精度（有效位数中的位数）。

标签： c floating-point

【解决方案1】：

有没有什么方法可以解释 IEEE 754 32 位浮点的精度舍入误差，它在不需要大量条件语句的情况下阻止 Out 达到 0？

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 是亚正常的，则只需返回 0.0f。

这类似于

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

【讨论】：

Subnormals 和 Out!=0 是红鲱鱼。标题询问到最终值，其他值也会出现同样的问题：Out 从In 得到 1 ULP，将其乘以 0.9（或 (.5, 1) 中的其他值）得到更多大于 In 的 ULP 的 ½，并将其添加到 In 会产生 Out。循环应该写得更通用，如float In = TargetValue; for (float Out = StartValue; Out != In;) Exp_Filt(In, &Out, 0.1f, 1.0f);
在循环（或重复调用的函数）内使用fesetround 可能会产生可怕的性能。在许多处理器上，更改浮点舍入模式会在一段时间内阻塞指令队列。
@EricPostpischilq 同意更改舍入模式有时间影响，但代码正在调用 expf()。我怀疑 OP 的时间问题不是这里的驱动问题。
@chux 一般来说，我喜欢你的回答，但为了性能，我会简化它，只计算(*Out - In) >= FLT_MIN，否则设置*Out = In;。
@Adam 我想你想要绝对值或fabsf(*Out - In) >= FLT_MIN。另外，请注意边缘情况。我推荐fabsf(*Out - In) > FLT_MIN。

【解决方案2】：

因为你正在失去精度，你应该以更高的精度进行中间计算（例如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) [once] 确实会收敛到零。

此外，根据给出的-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

编辑：原始帖子被 DV 化：

当你调用一个给定的函数时，如果函数的一个参数是float，它会被调用函数自动传递为double，并被解释为double 调用函数。

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

并且，在表达式中，float 在评估期间会自动提升为 double。

【讨论】：

(a) 默认参数提升（将float 转换为double）仅适用于函数声明中对应于... 的参数或没有原型的函数的参数。当一个参数在原型中声明为float 时，它就是float。 (b) 在评估期间float 值不会自动提升为double（尽管 C 实现可能会使用比标称类型更多的精度和范围进行评估）。 (c) 用这个答案中的 Exp_Filt 替换原来的 Exp_Filt 并不能解决问题中提出的问题。
@EricPostpischil 够公平的。但是，OP [担心] 失去精度 [和收敛]，因此将 args 更改为 double 会有所帮助。并且，该函数将 [稍微] 更快地运行 [因为 args 不必在表达式中提升/转换为 double。
不，OP 不关心丢失精度；这不是问题。而且，再一次，没有升级或转换为double。见上文 (b)。
关于编辑：问题中的函数不能解决问题。
C90 以后将 float 传递为 float 如果范围内有原型并且它没有 ... 省略号使问题复杂化。