一种基于AVR单片机的工频干扰滤除快速算法

liyf

一种基于AVR单片机的工频干扰滤除快速算法

摘要：提出了一种利用数字滤波器滤除工频干扰的快速算法。这种算法从AVR单片机内部硬件乘法器的特点出发，采用分配系统数法进行低通数字滤波器设计。经过VMLAB集成开发环境的仿真验证，算法速度快、代码效率高、滤波效果理想。     关键词：单片机 定点小数 FIR 工频干扰 分配系数法 VMLAB
工频干扰广泛存在各种工业现场中，其产生的途径主要包括输电馈线、照明设备、发动机以及各种电子仪器设备等。一般可以通过滤波电路消除工频干扰，但这必将增加硬件结构的复杂程度。实际上，还可以采用数字信号处理的相关算法，通过软件滤波器滤除工频干扰。软件滤波算法的采用，无疑会在简化电路结构的同时，使系统的硬件资源得到更加充分的利用，并达到降低产品成本的要求。
AVR单片机是Atmel公司生产的8位精简指令集（RISC）单片机。与同类单片机相比，在运算速度、外设资源、灵活性等方面性能均衡，性价比较高。AVR单片机适合C语言开发，Mega系列AVR单片机还有一个内部硬件乘法器单元。这些特点都为软件滤波器的实现提供了极大的便利。
1 滤波算法
常用的单片机滤波算法包括中值滤波、增色值滤波、滑动平均以及复合滤波算法等。工频干扰的频率范围在50Hz附近，可以采用一个截止频率远低于50Hz的低通滤波器来滤除工频干扰。

    假设输入信号x(t)，输出信号y(t)，则一个RC低通滤波器表达式为：τ[dy(t)]/dt+y(t)=x(t),其中τ为时间常数。
连续时间信号经过采样后成为离散时间信号，低通滤波器的表达式也变为：
τ[y(t)-y(k-1)]/ Δt+y(k)=x(t)     (1)
Δt为采样时间间隔，k为归一化时间。
由（1）式可得：
y(k)=[1/(1+τ/Δt)]x(k)+(τ/Δt)/(1+τ/Δt)[y(k-1)    (2)
令a=(τ/Δt)(1+τ/Δt)，代入（2）式，
y(k)=(1-a)·x(k)+a·y(k-1)     (3)
对（3）式进行z变换，可得：
Y(z)=(1-a) ·X(z)+a·z-1·Y(z)
所以系统的传递函数为：

假设采样频率Fs=500Hz,a分别取0.8、0.85、0.9、0.95，代入（4）式，利用matlab画出频率响应曲线，如图1。其中50Hz频率对应的幅度衰减见表1。

表1 50Hz幅度衰减

a	衰减
0.8	-9.4dB
0.85	-11.9dB
0.9	-15.5dB
0.95	-21.6dB

2 定点小数表达方式
为了精确构造数字滤波器，经常要用到浮点数据和系统。在进行浮点数乘法运算时，针对AVR单片机设计的C编译器例如AVR-GCC，需要加入额外的数学库函数进行编译。而这会使编译后程序的代码量增加、处理时间加长、处理器的开销也随之增加。为了更大限度地降低系统开销，提高程序效率，采用定点小数表示形式进行乘法运算是最佳选择。
AVR单片机是一种8位精简指令集（RISC）单片机。其中megaAVR系列内部都带有一个硬件乘法器，计算一次8位乘8位的定点乘法只需2个时钟周期。因此采用8位定点采样数据乘以8位系数的定点乘法方式完成滤波器算法是最高效的。

    低通滤波处理涉及的运算形式为一个纯小数系数和一个已知数据相乘再相加。因此将系数采用定点小数的表示形式，对于提高算法速度是至关重要的。
可以定义一种8位定点小数表示形成——Q8，其各位权系数如下：


Q8数的表示范围从0到1-2 -8=0.99609375，每两个数之间的间隔是2 -8（0.00390625），其所能表示的纯小数共有2 8=256个。例如11011000就表示2 -1+2 -2+2 -4+2 -5=0.84375，而11011001就是表示2 -1+2 -2+2 -4+2 -5+2 -8=0.84765625，因此0.84375和0.84765625之间的纯小数只能用这两个数中的一个近似表示了。这对于乘法计算的精度有一定的影响，但是由于滤波公式（3）中的系数a和(1-a)都是常数，在整体性能稳定的情况下，系数微小的不确定性对滤波器整体性能并没有太大的影响。
3 分配系数法原理
从（3）式可知，滤波算法可以用迭代计算实现，为保证每个新的输出值都可以作为下次计算的输入值，必须使输出值和输入值的位宽度一致。必须使输出值和输入值的位宽度一致。AVR单片机内部硬件乘法器的输出结果为16位，两次乘法运算的结果还要进行加法运算，其结果很有可能超过16位宽度。如果要进行迭代计算，就要将乘加运算的结果转化成8位表示方式。一种解决方法是用查表法实现乘法计算，这样运算结果就直接表达成8位定点数形式，不用进行表示方式的转化，但是这种方案要占用额外的硬件在座空间构造一张查找表。

    可以从逆向进行思考：由（3）式可知，每个新的输出值y(k)都与上一次的输出值y(k-1) 和新的输入值x(k)有关。y(k-1)和x(k)都是8位的，因此最大值为0xFF。为了使a×y(k-1)+(1-a) ×x(k)不超过0xFFFF，两个系统a和（1-a）的和不能超过0xFFFF/0xFF=0x101。实际上，a+(1-a)等于"1"，因此这里的0x101就可以看作“1”。如果取a=0.9，那么对应地将0x101平均分成10份，取其中的9份，即0x101×0.9近似等于0xE7，相应地0.1就等于0x101-0xE7=0x1A。这里的0xE7可以近似被认为是0.9的一种定点Q8数表示形式，而0.1的定点Q8数表示形式就是0x1A。由于滤波器系数a和（1-a）采用了Q8数的表示形式这种将16位乘加运算结果转化为8定点数表示形式的工作就变得各简单了，只需通过移位运算，取y(k)的高8位即可，对应的C语言代码为：
y(k)=（char）(y(k)>>8)
在C语言编程处理中，并不需要建立一个数组来存储y(k)的值，而只需定义两个unsigned char型的变量分别存储y(k-1)和x(k)。当乘加计算a×y(k-1)+(1-a)×x（k）完成后，将结果转化为8位定点数形式，再将其赋值给y(k-1)所对应的变量即可。因此采用迭代方式进行乘加运算后，整个运算过程只需要两个变量和两个常数参加即可。

    通过这种处理，y(k)就可以作为计算下一次输出值y(k+1)的一个已知量，并继续与Q8数形式的滤波器系数相乘，得到新的输出值。这种处理方式简化了乘加运算的完成过程，节省了系统硬件资源，并降低了处理器开销。
4 采样时间的控制
采用单片机进行数字信号处理，一种有效而准确的数据采集方式就是通过计数器中断服务程序（ISR）控制AD对输入信号进行精确采样。但是（图2）中断服务程序（ISR）的开销影响了AD采样时间间隔的精确度，同时如果中断服务程序（ISR）的开销过大，必然导致AD的最高采样频率的降低。因此，要想获得精确的采样频率，就必须在尽量减少中断服务程序开销的前提下，适当调整计数器中断的时间间隔。这可以通过调整OCR0的预置数来完成。
5 算法流程图
滤波算法是通过中断服务程序（ISR）来完成的，整个应用程序的主函数main（）主要负责初始化计数器中断，并处理其它应用。整个程序的流程图如图3所示。
本算法的C语言代码（附录A）经过AVR-GCC编译器的编译后，“.text”段只有310个字节，大大节省了单片机的flash空间。

6 基于VMLAB的滤波系统仿真实现
VMLAB的全称为：Visual Micro Lab。它针对AVR系列单片机和ST62系列单片机设计，是一个单片机的虚拟原型框架，可以提供给用户一个真正意义上的虚拟微控制器（MCU）设计实验室。它具有强大的多窗口、多文件的编辑器，微控制器的集成开发环境，拥有一系列的集成开发工具，图形界面的调试器，混合模式的模拟-数字电路仿真器，代码质量检测器等。基于MCU，它可以仿真出包括模拟元器件在内的更多外围设备，并具有交互式器件模拟仿真功能。
假设有用信号2V大小的直流信号，工频干扰是峰峰值为1V，频率为50Hz的正弦波，建立单片机AD的输入信号表示形式如下：
2+0.5 sin(2π×50×t)
VMLAB通过工程文件来管理和控制各种仿真信息、硬件连接以及显示I/O电压波形等。根据本算法的特点，采用Atmega16作为目标单片机，时钟选为8MHz，建立工程文件。恰当设置OCR0等存储，使计数器比较匹配中断的时间间隔约为2ms，这样AD的采样频率Fs近似认为等于500Hz。经过仿真，对比结果如表3。

表3 Fs=500Hz时仿真结果对比

a	DA输出纹波峰峰值	衰减幅度	DA输出均值	收敛时间
0.8	0.36V	-8.8dB	1.99V	约为100ms
0.85	0.28V	-11.1dB	1.99V	约为120ms
0.9	0.2V	-14dB	1.98V	约为145ms
0.95	0.1V	-20dB	1.95V	约为210ms

从表3可以看出：随着α的增大，算法收敛的时间变长，同时50Hz对应的衰减幅度增加，衰减的幅度值和理论推导基本一致。另外，当a=0.95时，DA输出的均值变小。这主要是进行循环迭代运算时，需要将16位的变量转化为8位表示形式所导致的。在有用信号失真较小的情况下，为使滤波器达到降低工频干扰的最佳效果，必须恰当选择a值。经过以上的仿真试验可以发现，当a=0.9时，衰减幅度、DA输出均值和算法收敛时间表现比较均衡，可以作为一般情况下的选择值。

    将VMLAB中虚拟示波器的显示数据导出到一个*.cvs文件中，用matlab读出这些数据，并画出不同a值对应的输出响应，如图4。从图4可以清晰看出不同a值下算法的性能变化的大致走向。
将AD的采样间隔设置为4ms，对应的采样频率Fs就变为250Hz，其它条件不变。通过VMLAB进行仿真，对比结果如表4、图5。

表4 Fs=250Hz时不同a值仿真结果对比

a	DA输出纹波峰峰值	衰减幅度	DA输出均值	收敛时间
0.8	0.20V	-14.0dB	1.98V	约为220ms
0.85	0.16V	-15.9dB	1.98V	约为240ms
0.9	0.14V	-17.1dB	1.97V	约为270ms
0.95	0.06V	-24.4dB	1.95V	大于500ms

对比Fs=500Hz的情况，随着采样频率Fs降低，50Hz频率的幅度衰减值会逐渐增加。这主要是因为随着采样频率降低，低通滤波器的截至频率fc也随之降低，相应的滤波器在50Hz处的衰减也就越来越低。根据奈奎斯特低通采样定理，当采样频率小于100Hz时，由于信号频谱混叠，滤波器对50Hz信号的滤波效果将会变差。如果只是对缓变信号进行采样，采样频率比100Hz稍大即可。但是随着采样频率的降低，滤波算法的收敛时间也会增加。因此必须在算法的滤波性能和收敛时间上进行折衷考虑。
本文提出的分配系数法设计数字滤波器，算法速度快、代码效率高、滤波效果理想，是一种实用的数字滤波器设计方法，体现了将算法嵌入到具体硬件的思想。另一方面，将定点小数的表示形式进行适当扩展，这个算法还可以用于10位或16位AD转换精度的应用场合。

李小路

谢谢分享！