一阶低通数字滤波器及其应用

一阶低通模拟滤波器

图中$u_i$ui​为输入电压，$u_o$uo​为输出电压。

基础方法推导公式

设流过电阻的电流为$i(t)$i(t)，因此我们有以下公式：

$q(t) = Cu_o(t)$q(t)=Cuo​(t)
$i(t) = d\frac{q(t)}{dt}$i(t)=ddtq(t)​
$u_i(t) = Ri(t)+u_o(t)$ui​(t)=Ri(t)+uo​(t)
式中$q(t)$q(t)为电容存储的电荷。
由以上三个式子得出：
$u_i(t) = RC\frac{du_o(t)}{dt}+u_o(t)$ui​(t)=RCdtduo​(t)​+uo​(t)

若输入电压为角频率为$\omega$ω的正弦波，即$u_i(t)=A\sin(\omega t+\varphi)$ui​(t)=Asin(ωt+φ)，带入上式得到

$RC\frac{du_o(t)}{dt}+u_o(t)=A\sin(\omega t+\varphi)$RCdtduo​(t)​+uo​(t)=Asin(ωt+φ)

上式是一个一阶微分方程，下面求其次方程的通解：

$RC\frac{du_o(t)}{dt}+u_o(t)=0$RCdtduo​(t)​+uo​(t)=0
$\Rightarrow \frac{1}{u_o(t)}du_o(t)=-\frac{1}{RC}dt$⇒uo​(t)1​duo​(t)=−RC1​dt
$\Rightarrow \ln\left|u_o(t)\right|=-\frac{t}{RC}+C_0$⇒ln∣uo​(t)∣=−RCt​+C0​
$\Rightarrow u_o(t)=C_1e^{-\frac{t}{RC}}$⇒uo​(t)=C1​e−RCt​

下面求解微分方程的特解：

很容易知道微分方程的一个特解形式为$u_o(t)=B\sin(\omega t+\theta)$uo​(t)=Bsin(ωt+θ)，代入方程式有
$B\omega RC\cos(\omega t+\theta)+B\sin(\omega t+\theta)=A\sin(\omega t+\varphi)$BωRCcos(ωt+θ)+Bsin(ωt+θ)=Asin(ωt+φ)
$\Rightarrow B\sqrt{(\omega RC)^2+1}\sin(\omega t+\theta+\alpha) = A\sin(\omega t+\varphi)$⇒B(ωRC)2+1​sin(ωt+θ+α)=Asin(ωt+φ)其中$\tan \alpha=\omega RC$tanα=ωRC
从上式可以得出：
$B=\frac{A}{\sqrt{(\omega RC)^2+1}}$B=(ωRC)2+1​A​、$\theta=\varphi-\alpha$θ=φ−α
故此，微分方程的特解为：
$u_o(t)=\frac{A}{\sqrt{(\omega RC)^2+1}}\sin(\omega t+\varphi-\alpha)$uo​(t)=(ωRC)2+1​A​sin(ωt+φ−α)，其中$\tan \alpha=\omega RC$tanα=ωRC

从而我们有微分方程的通解为：
$u_o(t)=C_1e^{-\frac{t}{RC}}+\frac{A}{\sqrt{(\omega RC)^2+1}}\sin(\omega t+\varphi-\alpha)$uo​(t)=C1​e−RCt​+(ωRC)2+1​A​sin(ωt+φ−α)，其中$\tan \alpha=\omega RC$tanα=ωRC

当时间趋向于无穷大时，上式中第一项趋于0。在研究滤波算法时，一般研究上式得第二项$\frac{A}{\sqrt{(\omega RC)^2+1}}\sin(\omega t+\varphi-\alpha)$(ωRC)2+1​A​sin(ωt+φ−α)。我们先看幅值变化，当角频率变大时，第二项的幅值$\frac{A}{\sqrt{(\omega RC)^2+1}}$(ωRC)2+1​A​变小，也就是说它具有通低频阻高频的效果。再来看看角度的变化，第二项的初始相位比输入延后了$\alpha$α角度，因此该电路的低通滤波具有延迟的效果。
截至频率$f_c$fc​的定义：滤波器的输出幅值变为输入幅值的$\frac{\sqrt{2}}{2}$22​​时的频率即为截至频率，即$\frac{A}{\sqrt{(\omega RC)^2+1}}=\frac{\sqrt{2}}{2}A$(ωRC)2+1​A​=22​​A，因此有截至频率：$f_c=\frac{\omega}{2\pi}=\frac{1}{2\pi RC}$fc​=2πω​=2πRC1​

相位法推导

由一阶低通滤波器电路可以得出：
$\dot{U}_o=\frac{1}{1+j\omega RC}\dot{U}_i=\frac{\dot{U}_i}{\sqrt{(\omega RC)^2+1}}\angle{-\alpha}$U˙o​=1+jωRC1​U˙i​=(ωRC)2+1​U˙i​​∠−α
其中$\omega$ω为输入正弦电压$u_i(t)$ui​(t)的角频率、$\tan\alpha=\omega RC$tanα=ωRC。设输入正弦电压$u_i(t)=A\sin(\omega t+\varphi)$ui​(t)=Asin(ωt+φ)，则输出也为正弦波，且表达式为：
$u_o(t)=\frac{A}{\sqrt{(\omega RC)^2+1}}\sin(\omega t+\varphi-\alpha)$uo​(t)=(ωRC)2+1​A​sin(ωt+φ−α)，此式与基础方式推导出的公式是一模一样。

一阶低通数字滤波器

滤波器公式推导

在模拟滤波器的基础方法中，有这样的式子：
$u_i(t) = RC\frac{du_o(t)}{dt}+u_o(t)$ui​(t)=RCdtduo​(t)​+uo​(t)
我们需要对上式进行离散化，假设采样频率为$f_s$fs​，则采样周期为$T=\frac{1}{f_s}$T=fs​1​。设$u_i(n)$ui​(n)、$u_o(n)$uo​(n)为当前采样时刻的值，$u_o(n-1)$uo​(n−1)为上一次采样时刻的值。则对上式离散化后有：
$u_i(n) = RCf_s[u_o(n)-u_o(n-1)]+u_o(n)$ui​(n)=RCfs​[uo​(n)−uo​(n−1)]+uo​(n)
整理后得到：
$u_o(n)=\frac{RCf_s}{RCf_s+1}u_o(n-1)+\frac{1}{RCf_s+1}u_i(n)$uo​(n)=RCfs​+1RCfs​​uo​(n−1)+RCfs​+11​ui​(n)
上式右边的两个系数之和为1，故此一阶低通滤波器公式也可写成：
$u_o(n)=au_i(n)+(1-a)u_o(n-1)$uo​(n)=aui​(n)+(1−a)uo​(n−1)，其中$0<a<1$0<a<1。

下面我们来研究$a$a与截止频率和延迟的关系。
根据上节及本节我们有公式:
$a=\frac{1}{RCf_s+1}$a=RCfs​+11​
$f_c=\frac{1}{2\pi RC}$fc​=2πRC1​
$\alpha=\arctan(\omega RC)$α=arctan(ωRC)
由此我们可以导出：
$a=\frac{2\pi f_c}{2\pi f_c+f_s}$a=2πfc​+fs​2πfc​​或
$f_c=\frac{af_s}{2\pi(1-a)}$fc​=2π(1−a)afs​​
如果输入的频率为$f$f，则输出的延迟时间为：
$T_{delay}=\frac{\alpha}{2\pi f}=\frac{\arctan(\frac{2\pi(1-a)f}{af_s})}{2\pi f}$Tdelay​=2πfα​=2πfarctan(afs​2π(1−a)f​)​

滤波器实例

数据准备

unsigned short DataIn[82] = 
    {391,420,442,465,466,472,478,485,491,495,
     498,488,480,479,480,482,481,480,476,468,
     451,444,442,441,441,437,433,429,427,410,
     399,393,389,388,388,384,379,369,360,351,
     341,331,327,323,320,318,317,315,311,305,
     298,288,274,265,260,255,252,250,249,247,
     246,245,244,241,236,229,222,216,210,204,
     200,197,193,190,186,181,178,177,175,173,
     171,169}

以上数据是通过5k采样频率得到的一个直流有刷电机停止时候的电流采样数据。我们需要根据以上数据来估计电机的位置。下图是根据采样数据画出来的。

由于电机在旋转过程中会有换相，因此上图的波形出现了一定的波动。我们目标是根据数据来估计上图波动的次数。

设计方案

如上图所示，输入数据x通过一阶低通滤波器后得到低频的数据，然后原始数据与低通数据相减后得到高频数据（从图中可以看出波动频率在300~800Hz），由于一阶低通滤波器有延时效应，故此得到如上图所示的高通滤波器。因此我们还需要计算一阶低通滤波器的延迟时间。
对于上图设置一阶低通滤波器的截止频率为110Hz。即$f_c=110$fc​=110，而$f_s=5000$fs​=5000。故而：
$a=\frac{2\pi f_c}{2\pi f_c+f_s}=0.1214$a=2πfc​+fs​2πfc​​=0.1214
由于我们是利用单片机计算的，故此可以取$a=\frac{1}{8}$a=81​。验算一下截止频率：
$f_c=\frac{af_s}{2\pi(1-a)}=113.7Hz$fc​=2π(1−a)afs​​=113.7Hz
下面来计算低频段的延迟时间：

频率(Hz)	延迟时间(ms)	延迟采样周期数
10	1.40	6.98
20	1.39	6.93
40	1.35	6.73
60	1.29	6.44
80	1.22	6.10
100	1.15	5.74
110	1.11	5.56

从表中可以看出该一阶低通滤波器对于低频段的延迟周期数大概为5~7。因此我们取6。
从而我们得到数字滤波器模型为：
$u(n)=\frac{1}{8}x(n)+\frac{7}{8}u(n-1)$u(n)=81​x(n)+87​u(n−1)
$y(n)=u(n)-x(n-6)$y(n)=u(n)−x(n−6)
C语言代码如下：

#define DELAY_CNT 6
/* startFlag:当input是第一个数据时候，调用该函数时，startFlag应该为1
 *         当input不是第一个数据时候，调用该函数时，startFlag应该为0
*/
unsigned short Filter(unsigned short input,unsigned char startFlag)
{
    static unsigned short temp[DELAY_CNT];
    static filterHis = 0;
    unsigned short filterTemp = 0;
    unsigned short output;
    static unsigned char saveCnt = 0;
    
    if(startFlag == 1)
    {
        int i;
        for(i=0; i<DELAY_CNT ; i++)
        {
            temp[i] = input;
            filterHis = input;
        }
    }
    else
    {
        filterTemp = input + filterHis * 7;
   	    temp[saveCnt++] = input;
        if(saveCnt == DELAY_CNT )
        {
            output = 1000 + (temp[0]<<3) - filterTemp;
            saveCnt = 0;
        }
        else
        {
            output = 2000 + (temp[saveCnt]<<3) - filterTemp;
        }
        filterHis = filterTemp>>3;
    }
    return output;
}

结果

得到的滤波后的数据如上图所示，单片机可以通过直接比大小的方式，就可以数清电机换相的次数。