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

一阶低通滤波器

结论

其基本原理基于以下公式：
$$output[n]=\alpha \ast input[n]+(1-\alpha )\ast 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)=\frac{1}{rcs+1}=\frac{1}{r\cdot jwc+1}=\frac{1}{\tau s+1}$$
在自控中称为一阶惯性环节
转成时域方程
$$X(s)=Y(s)/H(s)=Y(s)(\tau s+1)$$
$$x(t)=\tau y^{\prime }(t)+y(t)$$
将导数拆开(使用一阶后向差分法，对上面微分方程进行离散化)
$$x(t)=\tau (\frac{y(t)-y(t-T)}{T})+y(t)$$
整理成可递归迭代函数
$$y(t)=(1-\frac{T}{T+\tau })\cdot y(t-T)+\frac{T}{T+\tau }x(t)$$
令$a=\frac{T}{T+\tau }$ 则可得一般表达式
$$y(t)=(1-a)y(t-T)+ax(t)$$

2. 截止频率 和 采样频率 推导

从模电课本即可知，截止频率处幅值衰减-3db。

$$f_c=\frac{1}{2\pi RC}=\frac{1}{2\pi \tau }$$
将$\tau =\frac{1}{2\pi f_c}$带入滤波器系数$a$ 可得：
$$a=\frac{T}{T+\tau }=\frac{1}{1+\frac{f}{2\pi f_c}}$$
其中，$f=\frac{1}{T}$为采样频率。

实现

#include <stdio.h>
#include <stdlib.h>// 一阶低通滤波器函数
float lowPassFilter(float input, float prevOutput, float alpha) {return alpha * input + (1 - alpha) * prevOutput;
}int main() {float input = 10.0;  // 输入值float prevOutput = 5.0;  // 上一次的输出值float alpha = 0.2;  // 滤波系数for(int i=0; i<20; i++){float output = lowPassFilter(input, prevOutput, alpha);prevOutput = output;printf("filter current result: %f\n", output);}system("pause");return 0;
}

运行结果

二阶低通滤波器

实现1

#include <stdio.h>
// #include </lib/gcc/x86_64-linux-gnu/9/math.h>
#include <math.h>// 二阶低通滤波器参数
#define SAMPLING_FREQ 1000  // 采样频率
#define CUTOFF_FREQ 100  // 截止频率// 计算滤波器系数
void calculateFilterCoefficients(double *a, double *b) {double omega = 2 * M_PI * CUTOFF_FREQ / SAMPLING_FREQ;double alpha = sin(omega) / (2 * 0.707);double beta = cos(omega);double a0 = 1 + alpha;double a1 = -2 * beta;double a2 = 1 - alpha;double b0 = (1 - beta) / 2;double b1 = 1 - beta;double b2 = (1 - beta) / 2;*a = a0;*(a + 1) = a1;*(a + 2) = a2;*b = b0;*(b + 1) = b1;*(b + 2) = b2;
}// 二阶低通滤波函数
double lowPassFilter(double input, double *prevInputs, double *prevOutputs, double *a, double *b) {double output = *b * input + *b * prevInputs[0] + *b * prevInputs[1] - *a * prevOutputs[0] - *a * prevOutputs[1];prevInputs[1] = prevInputs[0];prevInputs[0] = input;prevOutputs[1] = prevOutputs[0];prevOutputs[0] = output;return output;
}int main() {double a[3], b[3];calculateFilterCoefficients(a, b);double prevInputs[2] = {0};double prevOutputs[2] = {0};double input = 10;  // 输入值，可根据实际情况修改double filteredOutput = lowPassFilter(input, prevInputs, prevOutputs, a, b);printf("滤波后的输出: %f\n", filteredOutput);return 0;
}/**
如果你在使用gcc编译含数学函数的 C 程序时，出现undefined reference to 'sin'、undefined reference to 'cos'等错误，一般是由于缺少库造成的。因为在 Ubuntu 系统中，gcc的数学函数（如sin、cos等）是定义在libm.so里面的，而数学库不在默认路径下。通过添加-lm选项，就可以告诉编译器到正确的库中查找这些函数。注意：在使用cmake进行编译时，需要添加命令target_link_libraries(your_target_name m)来链接数学库，其中your_target_name是你的目标名称。
*/

经我实际验证ubuntu20，的math库在如下路径
dpkg -l | grep math

实现2

#include <stdio.h>
#include <math.h>// 二阶低通滤波器系数
typedef struct {double a0, a1, a2, b1, b2;
} FilterCoefficients;// 计算二阶低通滤波器系数
void calculateFilterCoefficients(double cutoffFrequency, double samplingFrequency, FilterCoefficients *coefficients) {double omega = 2.0 * 3.14159 * cutoffFrequency / samplingFrequency;double cosOmega = cos(omega);double sinOmega = sin(omega);double alpha = sinOmega / (2.0 * 0.707);double a0 =  1 + alpha;double a1 = -2 * cosOmega;double a2 =  1 - alpha;double b1 = -2 * cosOmega;double b2 =  1 - alpha;coefficients->a0 = 1.0 / a0;coefficients->a1 = a1 / a0;coefficients->a2 = a2 / a0;coefficients->b1 = b1 / a0;coefficients->b2 = b2 / a0;
}// 二阶低通滤波器函数
void secondOrderLowPassFilter(double input[], double output[], int length, FilterCoefficients coefficients) {output[0] = input[0];output[1] = coefficients.a0 * input[1] + coefficients.a1 * input[0] + coefficients.b1 * output[0];for (int i = 2; i < length; i++) {output[i] = coefficients.a0 * input[i] + coefficients.a1 * input[i - 1] + coefficients.a2 * input[i - 2]- coefficients.b1 * output[i - 1] - coefficients.b2 * output[i - 2];}
}int main() {double cutoffFrequency = 10.0;  // 截止频率double samplingFrequency = 50.0;  // 采样频率FilterCoefficients coefficients;calculateFilterCoefficients(cutoffFrequency, samplingFrequency, &coefficients);double input[] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0};double output[10];int length = 10;secondOrderLowPassFilter(input, output, length, coefficients);for (int i = 0; i < length; i++) {printf("Output[%d] = %f\n", i, output[i]);}return 0;
}

实验结果：

推导1

二阶数字低通滤波器的时域上的开环传递函数为：
$$G(s)=\frac{w_c^2}{s^2+2\xi w_{c}s+w_c^2}$$
其中$w_c$是低通滤波器的截止频率，$\xi$为阻尼比。
采用双线性变换法 将 线性非时变系统 在 连续时域系统的传递函数 转换 成线性且平移不变滤波器 在 离散时域的传递函数。
双线性变换法：$S=\frac{2}{T}\frac{Z-1}{Z+1}$
可得
$$G(s)=\frac{w_c^2{T}_{sw}^2(z^2+2z+1)}{z^2(w_c^2{T}_{sw}^2+4\xi w_c{T}_{sw}+4)+z(2w_c^2{T}_{sw}^2-8)+(w_c^2{T}_{sw}^2-4\xi w_c{T}_{sw}+4)}$$

令系数为：
$$b_0=w_c^2$$
$$a_0=\frac{4}{T^2}+\frac{4}{T}\xi w_c+w_c^2$$
$$a_1=2w_c^2-\frac{8}{T^2}$$
$$a_2=w_2-2\xi w_c\frac{2}{T}+\frac{4}{T^2}$$
等效为：
$$b_0=w_c^{2}T$$
$$a_0=w_c^2{T}^2+4\xi w_{c}T+4$$
$$a_1=2w_c^2{T}^2-8$$
$$a_2=w_2{T}^2-4\xi w_{c}T+4$$

根据Z变换的时移性质
分母：$\frac{Y(n-1)}{Y(n)}=z^{-1}$, $\frac{Y(n-2)}{Y(n)}=z^{-2}$;
分子：$\frac{X(n-1)}{X(n)}=z^{-1}$，$\frac{X(n-2)}{X(n)}=z^{-2}$;
进一步可得：
$$Y(n)=\frac{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}$$

其中，$T_{sw}$为数据采样频率周期，$\xi$为阻尼比，取值0.707， $w_c=2\pi f_c$ (截止频率)

推导2

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