C/C++编程-算法学习-数字滤波器

打印 上一主题 下一主题

主题 840|帖子 840|积分 2520

一阶低通滤波器

结论

其基本原理基于以下公式:
                                         o                            u                            t                            p                            u                            t                            [                            n                            ]                            =                            α                            ∗                            i                            n                            p                            u                            t                            [                            n                            ]                            +                            (                            1                            −                            α                            )                            ∗                            o                            u                            t                            p                            u                            t                            [                            n                            −                            1                            ]                                  output[n] = α * input[n] + (1 - α) * output[n - 1]                     output[n]=α∗input[n]+(1−α)∗output[n−1]


  • output[n] 是当前的输出值
  • input[n] 是当前的输入值
  • α是滤波系数,取值范围在 0 到 1 之间。α值越大,对输入的相应越迅速,但滤波结果相对较弱;α值越小,滤波结果越强,但对输入的相应越慢
  • output[n - 1] 是上一次的输出值
    例如,如果 alpha = 0.1,输入值在短时间内快速变化,由于 (1 - alpha) 的权重较大,上一次的输出值对当前输出值的影响较大,从而起到平滑和抑制高频变化的作用。
推导1

1. 基本公式推导

对应电路模子一阶RC滤波器
                                         Y                            (                            s                            )                            /                            X                            (                            s                            )                            =                            H                            (                            s                            )                            =                                       1                                           r                                  c                                  s                                  +                                  1                                                 =                                       1                                           r                                  ⋅                                  j                                  w                                  c                                  +                                  1                                                 =                                       1                                           τ                                  s                                  +                                  1                                                       Y(s)/X(s)=H(s) = \frac{1}{rcs +1} = \frac{1}{r·jwc + 1} =\frac{1}{\tau s + 1}                     Y(s)/X(s)=H(s)=rcs+11​=r⋅jwc+11​=τs+11​
在自控中称为一阶惯性环节
转成时域方程
                                         X                            (                            s                            )                            =                            Y                            (                            s                            )                            /                            H                            (                            s                            )                            =                            Y                            (                            s                            )                            (                            τ                            s                            +                            1                            )                                  X(s) = Y(s)/H(s) = Y(s)(\tau s +1)                     X(s)=Y(s)/H(s)=Y(s)(τs+1)
                                         x                            (                            t                            )                            =                            τ                                       y                               ′                                      (                            t                            )                            +                            y                            (                            t                            )                                  x(t) = \tau y'(t) + y(t)                     x(t)=τy′(t)+y(t)
将导数拆开(使用一阶后向差分法,对上面微分方程进行离散化)
                                         x                            (                            t                            )                            =                            τ                            (                                                   y                                  (                                  t                                  )                                  −                                  y                                  (                                  t                                  −                                  T                                  )                                          T                                      )                            +                            y                            (                            t                            )                                  x(t) = \tau (\frac {y(t) - y(t-T)}{T}) + y(t)                     x(t)=τ(Ty(t)−y(t−T)​)+y(t)
整理成可递归迭代函数
                                         y                            (                            t                            )                            =                            (                            1                            −                                       T                                           T                                  +                                  τ                                                 )                            ⋅                            y                            (                            t                            −                            T                            )                            +                                       T                                           T                                  +                                  τ                                                 x                            (                            t                            )                                  y(t) = (1-\frac {T}{T+\tau})·y(t-T) + \frac{T}{T+\tau}x(t)                     y(t)=(1−T+τT​)⋅y(t−T)+T+τT​x(t)
令                                   a                         =                                   T                                       T                               +                               τ                                                 a = \frac{T}{T+\tau}                  a=T+τT​ 则可得一般表达式
                                         y                            (                            t                            )                            =                            (                            1                            −                            a                            )                            y                            (                            t                            −                            T                            )                            +                            a                            x                            (                            t                            )                                  y(t) = (1-a)y(t-T)+ax(t)                     y(t)=(1−a)y(t−T)+ax(t)
2. 截止频率 和 采样频率 推导

从模电课本即可知,截止频率处幅值衰减-3db。
                                                    f                               c                                      =                                       1                                           2                                  π                                  R                                  C                                                 =                                       1                                           2                                  π                                  τ                                                       f_c = \frac{1}{2\pi RC} = \frac{1}{2\pi \tau}                     fc​=2πRC1​=2πτ1​
将                                   τ                         =                                   1                                       2                               π                                           f                                  c                                                            \tau = \frac{1}{2\pi f_c}                  τ=2πfc​1​带入滤波器系数                                   a                              a                  a 可得:
                                         a                            =                                       T                                           T                                  +                                  τ                                                 =                                       1                                           1                                  +                                               f                                                   2                                        π                                                       f                                           c                                                                                              a = \frac{T}{T+ \tau} = \frac{1}{1+\frac{f}{2\pi f_c}}                     a=T+τT​=1+2πfc​f​1​
其中,                                   f                         =                                   1                            T                                       f = \frac{1}{T}                  f=T1​为采样频率。
实现

  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. // 一阶低通滤波器函数
  4. float lowPassFilter(float input, float prevOutput, float alpha) {
  5.     return alpha * input + (1 - alpha) * prevOutput;
  6. }
  7. int main() {
  8.     float input = 10.0;  // 输入值
  9.     float prevOutput = 5.0;  // 上一次的输出值
  10.     float alpha = 0.2;  // 滤波系数
  11.     for(int i=0; i<20; i++)
  12.     {
  13.         float output = lowPassFilter(input, prevOutput, alpha);
  14.         prevOutput = output;
  15.         printf("filter current result: %f\n", output);
  16.     }
  17.     system("pause");
  18.     return 0;
  19. }
复制代码
运行结果

二阶低通滤波器

实现1

  1. #include <stdio.h>
  2. // #include </lib/gcc/x86_64-linux-gnu/9/math.h>
  3. #include <math.h>
  4. // 二阶低通滤波器参数
  5. #define SAMPLING_FREQ 1000  // 采样频率
  6. #define CUTOFF_FREQ 100  // 截止频率
  7. // 计算滤波器系数
  8. void calculateFilterCoefficients(double *a, double *b) {
  9.     double omega = 2 * M_PI * CUTOFF_FREQ / SAMPLING_FREQ;
  10.     double alpha = sin(omega) / (2 * 0.707);
  11.     double beta = cos(omega);
  12.     double a0 = 1 + alpha;
  13.     double a1 = -2 * beta;
  14.     double a2 = 1 - alpha;
  15.     double b0 = (1 - beta) / 2;
  16.     double b1 = 1 - beta;
  17.     double b2 = (1 - beta) / 2;
  18.     *a = a0;
  19.     *(a + 1) = a1;
  20.     *(a + 2) = a2;
  21.     *b = b0;
  22.     *(b + 1) = b1;
  23.     *(b + 2) = b2;
  24. }
  25. // 二阶低通滤波函数
  26. double lowPassFilter(double input, double *prevInputs, double *prevOutputs, double *a, double *b) {
  27.     double output = *b * input + *b * prevInputs[0] + *b * prevInputs[1] - *a * prevOutputs[0] - *a * prevOutputs[1];
  28.     prevInputs[1] = prevInputs[0];
  29.     prevInputs[0] = input;
  30.     prevOutputs[1] = prevOutputs[0];
  31.     prevOutputs[0] = output;
  32.     return output;
  33. }
  34. int main() {
  35.     double a[3], b[3];
  36.     calculateFilterCoefficients(a, b);
  37.     double prevInputs[2] = {0};
  38.     double prevOutputs[2] = {0};
  39.     double input = 10;  // 输入值,可根据实际情况修改
  40.     double filteredOutput = lowPassFilter(input, prevInputs, prevOutputs, a, b);
  41.     printf("滤波后的输出: %f\n", filteredOutput);
  42.     return 0;
  43. }
  44. /**
  45. 如果你在使用gcc编译含数学函数的 C 程序时,出现undefined reference to 'sin'、undefined reference to 'cos'等错误,一般是由于缺少库造成的。因为在 Ubuntu 系统中,gcc的数学函数(如sin、cos等)是定义在libm.so里面的,而数学库不在默认路径下。通过添加-lm选项,就可以告诉编译器到正确的库中查找这些函数。
  46. 注意:在使用cmake进行编译时,需要添加命令target_link_libraries(your_target_name m)来链接数学库,其中your_target_name是你的目标名称。
  47. */
复制代码
经我现实验证ubuntu20,的math库在如下路径
dpkg -l | grep math



实现2

  1. #include <stdio.h>
  2. #include <math.h>
  3. // 二阶低通滤波器系数
  4. typedef struct {
  5.     double a0, a1, a2, b1, b2;
  6. } FilterCoefficients;
  7. // 计算二阶低通滤波器系数
  8. void calculateFilterCoefficients(double cutoffFrequency, double samplingFrequency, FilterCoefficients *coefficients) {
  9.     double omega = 2.0 * 3.14159 * cutoffFrequency / samplingFrequency;
  10.     double cosOmega = cos(omega);
  11.     double sinOmega = sin(omega);
  12.     double alpha = sinOmega / (2.0 * 0.707);
  13.     double a0 =  1 + alpha;
  14.     double a1 = -2 * cosOmega;
  15.     double a2 =  1 - alpha;
  16.     double b1 = -2 * cosOmega;
  17.     double b2 =  1 - alpha;
  18.     coefficients->a0 = 1.0 / a0;
  19.     coefficients->a1 = a1 / a0;
  20.     coefficients->a2 = a2 / a0;
  21.     coefficients->b1 = b1 / a0;
  22.     coefficients->b2 = b2 / a0;
  23. }
  24. // 二阶低通滤波器函数
  25. void secondOrderLowPassFilter(double input[], double output[], int length, FilterCoefficients coefficients) {
  26.     output[0] = input[0];
  27.     output[1] = coefficients.a0 * input[1] + coefficients.a1 * input[0] + coefficients.b1 * output[0];
  28.     for (int i = 2; i < length; i++) {
  29.         output[i] = coefficients.a0 * input[i] + coefficients.a1 * input[i - 1] + coefficients.a2 * input[i - 2]
  30.                     - coefficients.b1 * output[i - 1] - coefficients.b2 * output[i - 2];
  31.     }
  32. }
  33. int main() {
  34.     double cutoffFrequency = 10.0;  // 截止频率
  35.     double samplingFrequency = 50.0;  // 采样频率
  36.     FilterCoefficients coefficients;
  37.     calculateFilterCoefficients(cutoffFrequency, samplingFrequency, &coefficients);
  38.     double input[] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0};
  39.     double output[10];
  40.     int length = 10;
  41.     secondOrderLowPassFilter(input, output, length, coefficients);
  42.     for (int i = 0; i < length; i++) {
  43.         printf("Output[%d] = %f\n", i, output[i]);
  44.     }
  45.     return 0;
  46. }
复制代码
实验结果:

推导1

二阶数字低通滤波器的时域上的开环传递函数为:
                                         G                            (                            s                            )                            =                                                   w                                  c                                  2                                                                   s                                     2                                              +                                  2                                  ξ                                               w                                     c                                              s                                  +                                               w                                     c                                     2                                                                   G(s) = \frac{w_c^2}{s^2+2\xi w_c s + w_c^2}                     G(s)=s2+2ξwc​s+wc2​wc2​​
其中                                             w                            c                                       w_c                  wc​是低通滤波器的截止频率,                                   ξ                              \xi                  ξ为阻尼比。
采用双线性变换法 将 线性非时变体系 在 一连时域体系的传递函数 转换 成线性且平移稳定滤波器 在 离散时域的传递函数。
双线性变换法:                                   S                         =                                   2                            T                                                       Z                               −                               1                                                 Z                               +                               1                                                 S = \frac{2}{T}\frac{Z-1}{Z+1}                  S=T2​Z+1Z−1​
可得
                                         G                            (                            s                            )                            =                                                                w                                     c                                     2                                                           T                                                   s                                        w                                                  2                                              (                                               z                                     2                                              +                                  2                                  z                                  +                                  1                                  )                                                                   z                                     2                                              (                                               w                                     c                                     2                                                           T                                                   s                                        w                                                  2                                              +                                  4                                  ξ                                               w                                     c                                                           T                                                   s                                        w                                                           +                                  4                                  )                                  +                                  z                                  (                                  2                                               w                                     c                                     2                                                           T                                                   s                                        w                                                  2                                              −                                  8                                  )                                  +                                  (                                               w                                     c                                     2                                                           T                                                   s                                        w                                                  2                                              −                                  4                                  ξ                                               w                                     c                                                           T                                                   s                                        w                                                           +                                  4                                  )                                                       G(s) = \frac{w_c^2T_{sw}^2(z^2+2z+1)}{z^2(w_c^2T_{sw}^2 + 4\xi w_cT_{sw}+4)+z(2w_c^2T_{sw}^2-8)+(w_c^2T_{sw}^2-4\xi w_c T_{sw} + 4)}                     G(s)=z2(wc2​Tsw2​+4ξwc​Tsw​+4)+z(2wc2​Tsw2​−8)+(wc2​Tsw2​−4ξwc​Tsw​+4)wc2​Tsw2​(z2+2z+1)​
令系数为:
                                                    b                               0                                      =                                       w                               c                               2                                            b_0 = w_c^2                     b0​=wc2​
                                                    a                               0                                      =                                       4                                           T                                  2                                                 +                                       4                               T                                      ξ                                       w                               c                                      +                                       w                               c                               2                                            a_0 = \frac{4}{T^2 }+\frac{4}{T}\xi w_c +w_c^2                     a0​=T24​+T4​ξwc​+wc2​
                                                    a                               1                                      =                            2                                       w                               c                               2                                      −                                       8                                           T                                  2                                                       a_1 = 2w_c^2-\frac{8}{T^2}                     a1​=2wc2​−T28​
                                                    a                               2                                      =                                       w                               2                                      −                            2                            ξ                                       w                               c                                                 2                               T                                      +                                       4                                           T                                  2                                                       a_2 = w_2-2\xi w_c \frac{2}{T}+\frac{4}{T^2}                     a2​=w2​−2ξwc​T2​+T24​
等效为:
                                                    b                               0                                      =                                       w                               c                               2                                      T                                  b_0 = w_c^2T                     b0​=wc2​T
                                                    a                               0                                      =                                       w                               c                               2                                                 T                               2                                      +                            4                            ξ                                       w                               c                                      T                            +                            4                                  a_0 = w_c^2T^2 + 4\xi w_c T + 4                     a0​=wc2​T2+4ξwc​T+4
                                                    a                               1                                      =                            2                                       w                               c                               2                                                 T                               2                                      −                            8                                  a_1 = 2w_c^2T^2 - 8                     a1​=2wc2​T2−8
                                                    a                               2                                      =                                       w                               2                                                 T                               2                                      −                            4                            ξ                                       w                               c                                      T                            +                            4                                  a_2 = w_2T^2-4\xi w_c T+ 4                     a2​=w2​T2−4ξwc​T+4
根据Z变换的时移性子
分母:                                                        Y                               (                               n                               −                               1                               )                                                 Y                               (                               n                               )                                            =                                   z                                       −                               1                                                 \frac {Y(n-1)}{Y(n)} = z^{-1}                  Y(n)Y(n−1)​=z−1,                                                         Y                               (                               n                               −                               2                               )                                                 Y                               (                               n                               )                                            =                                   z                                       −                               2                                                 \frac{Y(n-2)}{Y(n)} = z^{-2}                  Y(n)Y(n−2)​=z−2;
分子:                                                        X                               (                               n                               −                               1                               )                                                 X                               (                               n                               )                                            =                                   z                                       −                               1                                                 \frac{X(n-1)}{X(n)} = z^{-1}                  X(n)X(n−1)​=z−1,                                                        X                               (                               n                               −                               2                               )                                                 X                               (                               n                               )                                            =                                   z                                       −                               2                                                 \frac{X(n-2)}{X(n)} = z^{-2}                  X(n)X(n−2)​=z−2;
进一步可得:
                                         Y                            (                            n                            )                            =                                                                b                                     0                                              X                                  (                                  n                                  )                                  +                                               b                                     1                                              X                                  (                                  n                                  −                                  1                                  )                                  +                                               b                                     2                                              X                                  (                                  n                                  −                                  2                                  )                                  −                                               a                                     1                                              Y                                  (                                  n                                  −                                  1                                  )                                  −                                               a                                     2                                              Y                                  (                                  n                                  −                                  2                                  )                                                      a                                  0                                                       Y(n) = \frac {b_0X(n) + b_1X(n-1) + b_2X(n-2)-a_1Y(n-1)-a_2Y(n-2)}{a_0}                     Y(n)=a0​b0​X(n)+b1​X(n−1)+b2​X(n−2)−a1​Y(n−1)−a2​Y(n−2)​
其中,                                             T                                       s                               w                                                 T_{sw}                  Tsw​为数据采样频率周期,                                   ξ                              \xi                  ξ为阻尼比,取值0.707,                                              w                            c                                  =                         2                         π                                   f                            c                                       w_c = 2\pi f_c                  wc​=2πfc​ (截止频率)
推导2

【推导方法1】
【推导方法2】

免责声明:如果侵犯了您的权益,请联系站长,我们会及时删除侵权内容,谢谢合作!更多信息从访问主页:qidao123.com:ToB企服之家,中国第一个企服评测及商务社交产业平台。

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有账号?立即注册

x
回复

使用道具 举报

0 个回复

倒序浏览

快速回复

您需要登录后才可以回帖 登录 or 立即注册

本版积分规则

鼠扑

金牌会员
这个人很懒什么都没写!

标签云

快速回复 返回顶部 返回列表