一阶低通滤波的C语言实现

  • 0 引言
  • 1 一阶低通滤波器模型
  • 2 matlab 实现
    • 2.1 matlab 代码
    • 2.2 总结
  • 3 c语言实现
  • 4 matlab 查看波形频率（快速傅里叶变换,FFT）

0 引言

一阶低通滤波器（Low Pass Filter，LPF），核心参数为截止频率fc，该算法可以保留截止频率以内的信号，而衰减截止频率之外的信号。主要用于去除高频噪声。

1 一阶低通滤波器模型

一阶低通滤波公式如下：

也可以写作：

其中：

参数说明：y(n)为本次滤波输出值，y(n-1)为上次滤波输出值，x(n)为本次采样值。Ts为采样周期，fc为截止频率。α范围为[0,1]

2 matlab 实现

我们假设，现在有一个信号，它包含了频率为1Hz(幅值为3)和4Hz(幅值为1)的两个正弦波。

如图所示：黄色为叠加后的信号。
我们假定采样周期为0.02s,且需要保留1Hz的信号(那么截止频率应该在1Hz~4Hz之间，先取1Hz)。

则

即

2.1 matlab 代码

t=(0:0.02:5);x=3*sin(2*pi*1*t)+1*sin(2*pi*4*t);%采集的混合信号plot(t,x);hold on fc=2;%截止频率Ts=0.02;%采样周期b=2*pi*fc*Ts;alpha=b/(b+1);%alphay=zeros(1,251);y(1)=x(1)for n= 2:251y(n)=y(n-1)+alpha*(x(n)-y(n-1));%一阶滤波endplot(t,y);

我们将滤波后的波形和理论的1Hz的波形对比：

你会发现：

  1. 确实有效果
  2. 波形幅度变小了
  3. 波形滞后了

因此，一阶低通滤波，存在滞后的现象。那么我们可以适当修改一下α，比如，将截止频率设置为2Hz，再看看效果。

你会发现：

  • 幅值变化不大
  • 相位滞后不大

2.2 总结

对于一阶低通滤波器，截止频率应在有效信号频率和杂波频率之间。
越接近有效信号频率，则滤波效果越强，但存在相位更加滞后，幅度衰减越厉害的情况。
因此，应该多次设值，从而寻找一个能够兼顾响应速度与滤波效果的值。

3 c语言实现

虽然在网上看到别人写的c语言实现，但基本是原理说明，搬过来没办法直接使用，还要修修改改的，就很难受。
我这里模拟采集信号的过程，写一个c语言的实现。主要添加了一个模拟输入的效果,数据其实就是从matlab里面复制出来的：x=3*sin(2*pi*1*t)+1*sin(2*pi*4*t)。

#include #include /*************************** 模拟输入数据 ********************************/float input_data[251] = {0.00, 0.86, 1.59, 2.10, 2.35, 2.35, 2.18, 1.94, 1.76, 1.73, 1.90, 2.26, 2.75,3.24, 3.63, 3.80, 3.70, 3.30, 2.68, 1.93, 1.18, 0.54, 0.11, -0.10, -0.11,0.00, 0.11, 0.10, -0.11, -0.54, -1.18, -1.93, -2.68, -3.30, -3.70, -3.80, -3.63, -3.24, -2.75, -2.26, -1.90, -1.73, -1.76, -1.94, -2.18, -2.35, -2.35,-2.10, -1.59, -0.86, 0.00, 0.86, 1.59, 2.10, 2.35, 2.35, 2.18, 1.94, 1.76, 1.73, 1.90, 2.26, 2.75, 3.24, 3.63, 3.80, 3.70, 3.30, 2.68, 1.93, 1.18,0.54, 0.11, -0.10, -0.11, 0.00, 0.11, 0.10, -0.11, -0.54, -1.18, -1.93, -2.68, -3.30, -3.70, -3.80, -3.63, -3.24, -2.75, -2.26, -1.90, -1.73,-1.76, -1.94, -2.18, -2.35, -2.35, -2.10, -1.59, -0.86, 0.00, 0.86, 1.59,2.10, 2.35, 2.35, 2.18, 1.94, 1.76, 1.73, 1.90, 2.26, 2.75, 3.24, 3.63, 3.80, 3.70, 3.30, 2.68, 1.93, 1.18, 0.54, 0.11, -0.10, -0.11, 0.00, 0.11,0.10, -0.11, -0.54, -1.18, -1.93, -2.68, -3.30, -3.70, -3.80, -3.63, -3.24,-2.75, -2.26, -1.90, -1.73, -1.76, -1.94, -2.18, -2.35, -2.35, -2.10, -1.59,-0.86, 0.00, 0.86, 1.59, 2.10, 2.35, 2.35, 2.18, 1.94, 1.76, 1.73, 1.90, 2.26,2.75, 3.24, 3.63, 3.80, 3.70, 3.30, 2.68, 1.93, 1.18, 0.54, 0.11, -0.10, -0.11,0.00, 0.11, 0.10, -0.11, -0.54, -1.18, -1.93, -2.68, -3.30, -3.70, -3.80, -3.63, -3.24, -2.75, -2.26, -1.90, -1.73, -1.76, -1.94, -2.18, -2.35, -2.35, -2.10, -1.59,-0.86, 0.00, 0.86, 1.59, 2.10, 2.35, 2.35, 2.18, 1.94, 1.76, 1.73, 1.90, 2.26, 2.75, 3.24, 3.63, 3.80, 3.70, 3.30, 2.68, 1.93, 1.18, 0.54, 0.11, -0.10, -0.11, 0.00, 0.11,0.10, -0.11, -0.54, -1.18, -1.93, -2.68, -3.30, -3.70, -3.80, -3.63, -3.24, -2.75,-2.26, -1.90, -1.73, -1.76, -1.94, -2.18, -2.35, -2.35, -2.10, -1.59, -0.86, 0.00};float get_data(void){static int i = 0;return (input_data[i++]);if (i == 251) //轮回i = 0;}float fc = 2.0f; //截止频率float Ts = 0.02f;//采样周期float pi = 3.14159f; //πfloat alpha = 0; //滤波系数/************************ 滤波器初始化 alpha *****************************/void low_pass_filter_init(void){float b = 2.0 * pi * fc * Ts;alpha = b / (b + 1);}float low_pass_filter(float value){static float out_last = 0; //上一次滤波值float out;/***************** 如果第一次进入，则给 out_last 赋值 ******************/static char fisrt_flag = 1;if (fisrt_flag == 1){fisrt_flag = 0;out_last = value;}/*************************** 一阶滤波 *********************************/out = out_last + alpha * (value - out_last);out_last = out;return out;}void main(void){float result[251];low_pass_filter_init();for (int i = 0; i < 251; i++){result[i] = low_pass_filter(get_data());printf("%f,", result[i]);}}

我这里使用的是VSCode，将终端输出的数据直接复制到matlab加点代码就可以观察了。（excel也行，但插入数据后会坨在一格里，要点击分列）。
效果如下：

4 matlab 查看波形频率（快速傅里叶变换,FFT）

有时候，我们想先分析我们的数据的频率组成，或者观察滤波效果，我们可以观察他们的频率。这里用的是快速傅里叶变换。
这里分析X=3*sin(2*pi*1*t)+1*sin(2*pi*4*t)，代码如下：

t=(0:0.02:5);% 时间向量X=3*sin(2*pi*1*t)+1*sin(2*pi*4*t);%采集的混合信号Ts = 0.02;% 采样周期Fs = 1/Ts;% 采样频率 L = length(t);% 获取信号长度n = 2^nextpow2(L);%首先从原始信号长度确定下一个 2 次幂的新输入长度。这将用尾随零填充信号 X 以改善 fft 的性能。%%%%%%%%%%%%%%%%%%%%%%% 快速傅里叶变换 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Y = fft(X,n);%信号的傅里叶变换P2 = abs(Y/n);%计算双侧频谱 P2P1 = P2(1:n/2+1);%基于 P2 和偶数信号长度 L 计算单侧频谱 P1。P1(2:end-1) = 2*P1(2:end-1);%双侧边变成单侧侧边，幅度叠加f = Fs*(0:(n/2))/n;%定义频域 f plot(f,P1) title('Single-Sided Amplitude Spectrum of X(t)')xlabel('f (Hz)')ylabel('|P1(f)|')

我们将之前代码的滤波输出（fc=2）直接赋值给X,进行观察：

可以看到，频率为4Hz的信号幅度衰减超过一半，而1Hz的衰减非常少。说明效果还是不错的。