一阶低通滤波的C语言实现(简单易移植)

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);%alpha

y=zeros(1,251);

y(1)=x(1)
for n= 2:251
    y(n)=y(n-1)+alpha*(x(n)-y(n-1));%一阶滤波
end

plot(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 <stdio.h>
#include <string.h>

/*************************** 模拟输入数据 ********************************/
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);%计算双侧频谱 P2

P1 = 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的衰减非常少。说明效果还是不错的。