快速和内存高效的移动平均计算
正如其他人所提到的,您应该考虑使用 IIR(无限脉冲响应)滤波器,而不是您现在使用的 FIR(有限脉冲响应)滤波器。还有更多,但乍一看,FIR 滤波器是作为显式卷积和带有方程的 IIR 滤波器实现的。
我在微控制器中经常使用的特定 IIR 滤波器是单极低通滤波器。这是一个简单的 RC 模拟滤波器的数字等效物。对于大多数应用程序,这些将比您使用的盒式过滤器具有更好的特性。我遇到的大多数箱式滤波器的使用是由于有人没有关注数字信号处理课程,而不是因为需要它们的特定特性。如果您只想衰减已知为噪声的高频,则单极低通滤波器会更好。在微控制器中实现数字化的最佳方法通常是:
过滤 <-- 过滤 + FF(新 - 过滤)
FILT 是一种持久状态。这是计算此过滤器所需的唯一持久变量。NEW 是过滤器使用此迭代更新的新值。FF 是过滤器分数,它调整过滤器的“重量”。看看这个算法,发现对于 FF = 0,滤波器是无限重的,因为输出永远不会改变。对于 FF = 1,它实际上根本没有过滤器,因为输出只是跟随输入。有用的值介于两者之间。在小型系统上,您选择 FF 为 1/2 N因此乘以 FF 可以通过右移 N 位来完成。例如,FF 可能是 1/16,乘以 FF 因此右移 4 位。否则这个过滤器只需要一个减法和一个加法,尽管数字通常需要比输入值更宽(在下面的单独部分中更多关于数值精度)。
我通常会比所需的速度更快地读取 A/D 读数,并级联应用这些过滤器中的两个。这是两个串联 RC 滤波器的数字等效值,在滚降频率之上衰减 12 dB/倍频程。但是,对于 A/D 读数,通常更相关的是通过考虑其阶跃响应来查看时域中的滤波器。这告诉您当您正在测量的事物发生变化时,您的系统会以多快的速度看到变化。
为了便于设计这些过滤器(这只意味着选择 FF 并决定要级联的数量),我使用了我的程序 FILTBITS。您可以为级联系列滤波器中的每个 FF 指定移位位数,它会计算阶跃响应和其他值。实际上,我通常通过我的包装脚本 PLOTFILT 运行它。这会运行 FILTBITS,它会生成一个 CSV 文件,然后绘制 CSV 文件。例如,这里是“PLOTFILT 4 4”的结果:
PLOTFILT 的两个参数意味着将有两个上述类型的过滤器级联。4 的值表示实现乘以 FF 的移位位数。因此,在这种情况下,两个 FF 值是 1/16。
红色迹线是单位阶跃响应,是要查看的主要内容。例如,这告诉您如果输入瞬时改变,组合滤波器的输出将在 60 次迭代中稳定到新值的 90%。如果您关心 95% 的稳定时间,那么您必须等待大约 73 次迭代,而 50% 的稳定时间只需等待 26 次迭代。
绿色迹线显示了单个全幅尖峰的输出。这使您对随机噪声抑制有所了解。看起来没有一个样本会导致输出的变化超过 2.5%。
蓝色迹线是对这个滤波器对白噪声的作用的主观感受。这不是一个严格的测试,因为无法保证在这次 PLOTFILT 运行中选择作为白噪声输入的随机数的确切内容是什么。只是给你一个粗略的感觉,它会被压扁多少,它有多光滑。
PLOTFILT,也许是 FILTBITS,以及许多其他有用的东西,特别是对于 PIC 固件开发,可在我的软件下载页面的 PIC 开发工具软件版本中找到。
添加了关于数值精度
我从评论中看到,现在有一个新的答案,即有兴趣讨论实现此过滤器所需的位数。请注意,乘以 FF 将在二进制点下方创建 Log 2 (FF) 新位。在小型系统上,FF 通常选择为 1/2 N,因此这种乘法实际上是通过右移 N 位来实现的。
因此,FILT 通常是一个定点整数。请注意,从处理器的角度来看,这不会改变任何数学运算。例如,如果您要过滤 10 位 A/D 读数并且 N = 4 (FF = 1/16),那么您需要在 10 位整数 A/D 读数以下 4 个小数位。大多数处理器,由于 10 位 A/D 读数,您将执行 16 位整数运算。在这种情况下,您仍然可以执行完全相同的 16 位整数运算,但从 A/D 读数左移 4 位开始。处理器不知道区别,也不需要。无论您认为它们是 12.4 定点还是真正的 16 位整数(16.0 定点),对整个 16 位整数进行数学运算都有效。
通常,如果您不想因数字表示而添加噪声,则需要在每个滤波器极点添加 N 位。在上面的示例中,第二个过滤器必须有 10+4+4 = 18 位才能不丢失信息。实际上,在 8 位机器上,这意味着您将使用 24 位值。从技术上讲,只有两个中的第二个极点需要更宽的值,但为了简化固件,我通常对滤波器的所有极点使用相同的表示,因此使用相同的代码。
通常我会编写一个子程序或宏来执行一个滤波器极点操作,然后将其应用于每个极点。是子程序还是宏取决于循环或程序存储器在该特定项目中是否更重要。无论哪种方式,我都使用一些暂存状态将 NEW 传递到子例程/宏中,它会更新 FILT,但也会将其加载到 NEW 所在的相同暂存状态。这使得应用多个极点变得容易,因为一个极点的更新 FILT 是下一个的新。当一个子程序时,有一个指向 FILT 的指针很有用,它在输出时更新为 FILT 之后。这样,如果多次调用,子程序会自动对内存中的连续过滤器进行操作。使用宏不需要指针,因为您传入地址以对每次迭代进行操作。
代码示例
以下是上述 PIC 18 的宏示例:
///////////////////////////////////////// ////////////////////////// // // 宏 FILTER 过滤 // // 用 NEWVAL 中的新值更新一个滤波器极点。NEWVAL 更新为 // 包含新的过滤值。 // // FILT 是过滤器状态变量的名称。假设为 24 位 // 广泛且在本地银行中。 // // 更新过滤器的公式为: // // 过滤 <-- 过滤 + FF(NEWVAL - 过滤) // // 乘以 FF 是通过 FILTBITS 位的右移来完成的。 // /宏过滤器 /写 dbankif lbankadr movf [arg 1]+0, w ;NEWVAL <-- NEWVAL - 过滤 subwf newval+0 movf [arg 1]+1, w subwfb newval+1 movf [arg 1]+2, w subwfb newval+2 /写 /loop n filtbits ; 每一位将 NEWVAL 右移一次 rlcf newval+2, w ;将 NEWVAL 右移一位 rrcf newval+2 rrcf newval+1 rrcf newval+0 /endloop /写 movf newval+0, w ; 将移位值添加到过滤器中并保存在 NEWVAL addwf [arg 1]+0, w movwf [参数 1]+0 movwf newval+0 movf newval+1, w addwfc [arg 1]+1, w movwf [参数 1]+1 movwf newval+1 movf newval+2, w addwfc [arg 1]+2, w movwf [参数 1]+2 movwf newval+2 /endmac
这是 PIC 24 或 dsPIC 30 或 33 的类似宏:
///////////////////////////////////////// ////////////////////////// // // 宏过滤器 ffbits // // 更新一个低通滤波器的状态。新的输入值在 W1:W0 // 并且要更新的过滤器状态由 W2 指向。 // // 更新后的过滤器值也会在 W1:W0 中返回,W2 将指向 // 到过滤器状态后的第一个内存。因此这个宏可以是 // 连续调用以更新一系列级联低通滤波器。 // // 过滤公式为: // // 过滤 <-- 过滤 + FF(新 - 过滤) // // 其中乘以 FF 是通过算术右移执行的 // FFBITS。 // // 警告:W3 已被丢弃。 // /宏过滤器 /var new ffbits integer = [arg 1] ;获取要移位的位数 /写 /write " ; 执行一极低通滤波,移位位 = " ffbits /写 ” ;” sub w0, [w2++], w0 ;新 - 过滤 --> W1:W0 subb w1, [w2--], w1 lsr w0, #[v ffbits], w0 ;将 W1:W0 中的结果右移 sl w1, #[- 16 ffbits], w3 ior w0, w3, w0 asr w1, #[v ffbits], w1 add w0, [w2++], w0 ;添加 FILT 以在 W1:W0 中生成最终结果 addc w1, [w2--], w1 mov w0, [w2++] ;将结果写入过滤器状态,提前指针 mov w1, [w2++] /写 /endmac
这两个示例都是使用我的PIC 汇编器预处理器作为宏实现的,它比任何一个内置的宏工具都更有能力。
如果您可以忍受两个项目平均数(即 2、4、8、16、32 等)的幂的限制,那么可以在没有专用除法的低性能微机上轻松有效地完成除法,因为它可以作为位移来完成。每个右移是 2 的一次幂,例如:
avg = sum >> 2; //divide by 2^2 (4)
要么
avg = sum >> 3; //divide by 2^3 (8)
等等。
使用 Olin Lathrop 已经在数字信号处理堆栈交换中描述的一阶 IIR 滤波器背后的数学进行了一些深入的分析(包括许多漂亮的图片。)这个 IIR 滤波器的等式是:
y[n]=αx[n]+(1−α)y[n−1]
这可以使用以下代码仅使用整数而不是除法来实现(可能需要一些调试,因为我是从内存中输入的。)
/**
* @details Implement a first order IIR filter to approximate a K sample
* moving average. This function implements the equation:
*
* y[n] = alpha * x[n] + (1 - alpha) * y[n-1]
*
* @param *filter - a Signed 15.16 fixed-point value.
* @param sample - the 16-bit value of the current sample.
*/
#define BITS 2 ///< This is roughly = log2( 1 / alpha )
short IIR_Filter(long *filter, short sample)
{
long local_sample = sample << 16;
*filter += (local_sample - *filter) >> BITS;
return (short)((*filter+0x8000) >> 16); ///< Round by adding .5 and truncating.
}
该过滤器通过将 alpha 的值设置为 1/K 来近似最后 K 个样本的移动平均值。在前面的代码中通过对 LOG2(K) 执行此#define
操作BITS
,即 K = 16 设置BITS
为 4,K = 4 设置BITS
为 2,等等。
(我将在收到更改后立即验证此处列出的代码,并在需要时编辑此答案。)
如果您不介意下采样,那么内存需求较少的真正移动平均滤波器(又名“boxcar 滤波器”)有一个答案。它被称为级联积分梳状滤波器(CIC)。这个想法是您有一个积分器,您可以在一段时间内获取差异,而节省内存的关键设备是通过下采样,您不必存储积分器的每个值。它可以使用以下伪代码来实现:
function out = filterInput(in)
{
const int decimationFactor = /* 2 or 4 or 8 or whatever */;
const int statesize = /* whatever */
static int integrator = 0;
static int downsample_count = 0;
static int ringbuffer[statesize];
// don't forget to initialize the ringbuffer somehow
static int ringbuffer_ptr = 0;
static int outstate = 0;
integrator += in;
if (++downsample_count >= decimationFactor)
{
int oldintegrator = ringbuffer[ringbuffer_ptr];
ringbuffer[ringbuffer_ptr] = integrator;
ringbuffer_ptr = (ringbuffer_ptr + 1) % statesize;
outstate = (integrator - oldintegrator) / (statesize * decimationFactor);
}
return outstate;
}
您的有效移动平均长度是decimationFactor*statesize
,但您只需要保留statesize
样本。statesize
显然,如果您的和是 2 的幂,那么您可以获得更好的性能decimationFactor
,因此除法和余数运算符被移位和掩码与取代。
后记:我同意 Olin 的观点,即您应该始终在移动平均滤波器之前考虑简单的 IIR 滤波器。如果您不需要 boxcar 滤波器的频率零点,则 1 极或 2 极低通滤波器可能工作得很好。
另一方面,如果您出于抽取目的进行过滤(采用高采样率输入并将其平均以供低速率过程使用),那么 CIC 滤波器可能正是您正在寻找的。(特别是如果您可以使用 statessize=1 并且只使用一个先前的积分器值就可以完全避免环形缓冲区)