文章目录
- 1、一阶低通滤波
- 2、一阶高通滤波
- 3、二阶低通滤波器
- 3.1 二阶RC低通滤波器的连续域数学模型
- 3.2 二阶RC低通滤波器的算法推导
- 3.3 matlab仿真
- 4、二阶高通滤波器
- 4.1 二阶RC高通滤波器的连续域数学模型
- 4.2 二阶RC高通滤波器的算法推导
- 4.3 matlab仿真
- 5、陷波滤波
- 6、带通滤波器
- 6.1 带通滤波器的连续域数学模型
- 6.2 带通滤波器的算法推导
- 6.3 matlab仿真
- 参考
写在前面:本文是对《RC滤波器数学公式推导及软件算法实现》的修正和完善,修正主要是针对二阶滤波的修正,并对所有滤波算法进行了matlab算法验证,利用FFT对滤波前和滤波后的信号进行了分析。
1、一阶低通滤波
差分离散表达式:
V
o
(
k
)
=
A
V
i
(
k
)
+
(
1
−
A
)
V
o
(
k
−
1
)
V_o(k)=AV_i(k)+(1-A)V_o(k-1)
Vo(k)=AVi(k)+(1−A)Vo(k−1)
其中
T
s
T_s
Ts为采样频率,
R
C
=
1
/
2
π
f
c
RC=1/2\pi f_c
RC=1/2πfc,于是已知截至频率,就可以得到系数,也可以由系数求得截至频率。
A
=
w
c
T
s
1
+
w
c
T
s
A=\frac{w_c T_s}{1+w_c T_s}
A=1+wcTswcTs
上式就是已知截至频率
f
c
f_c
fc,可得系数
A
A
A
matlab验证代码:
% 参数设置
fs = 1000; % 采样频率 1000 Hz
wc = 2 * pi * 10; % 截止频率 10 Hz
Ts = 1/fs; % 采样周期
A = (wc * Ts) / (1 + wc * Ts); % 计算系数 A
% 仿真输入信号 (正弦和噪声信号的叠加)
t = 0:Ts:1-Ts; % 离散时间向量
input_signal = sin(2*pi*5*t) + sin(2*pi*50*t) + 0.5*randn(size(t)); % 混合信号: 5 Hz + 50 Hz + 噪声
% 初始化输出信号
output_signal = zeros(size(input_signal));
% 初始条件设置
output_signal(1) = A * input_signal(1); % 假设初始输出为 0
% 递推公式计算一阶低通滤波器输出
for k = 2:length(input_signal)
output_signal(k) = A * input_signal(k) + (1 - A) * output_signal(k-1);
end
% FFT 分析
n = length(t); % 数据点数
f = (0:n-1)*(fs/n); % 频率向量
input_fft = abs(fft(input_signal)); % 输入信号的 FFT
output_fft = abs(fft(output_signal)); % 输出信号的 FFT
% 画图 - 输入输出信号的时域响应
figure;
subplot(2, 1, 1);
plot(t, input_signal);
title('输入信号 (时域)');
xlabel('时间 (s)');
ylabel('振幅');
subplot(2, 1, 2);
plot(t, output_signal);
title('滤波后信号 (时域)');
xlabel('时间 (s)');
ylabel('振幅');
% 画图 - 输入输出信号的频域响应
figure;
subplot(2, 1, 1);
plot(f(1:n/2), input_fft(1:n/2)); % 只显示正频率部分
title('输入信号 (频域)');
xlabel('频率 (Hz)');
ylabel('幅度');
subplot(2, 1, 2);
plot(f(1:n/2), output_fft(1:n/2)); % 只显示正频率部分
title('滤波后信号 (频域)');
xlabel('频率 (Hz)');
ylabel('幅度');
显示滤波前后信号的振幅随时间的变化,滤波后信号的高频噪声明显被抑制 。
FFT 分析显示滤波前信号的频谱具有 5 Hz 和 50 Hz 的显著频率成分,而经过低通滤波后,高频(50 Hz)分量被显著衰减,只剩下低频(5 Hz)的成分
2、一阶高通滤波
差分离散表达式:
V
o
(
k
)
=
A
V
0
(
k
−
1
)
+
A
(
V
i
(
k
)
−
V
i
(
k
−
1
)
)
V_o(k)=AV_0(k-1)+A(V_i(k)-V_i(k-1))
Vo(k)=AV0(k−1)+A(Vi(k)−Vi(k−1))
其中
T
s
T_s
Ts为采样频率,
R
C
=
1
/
2
π
f
c
RC=1/2\pi f_c
RC=1/2πfc,于是已知截至频率,就可以得到系数,也可以由系数求得截至频率。
A
=
1
1
+
w
c
T
s
A=\frac{1}{1+w_c T_s}
A=1+wcTs1
上式就是已知截至频率
f
c
f_c
fc,可得系数
A
A
A
matlab算法验证:
% 参数设置
fs = 1000; % 采样频率 1000 Hz
wc = 2 * pi * 50; % 截止频率 50 Hz
Ts = 1/fs; % 采样周期
A = 1 / (1 + wc * Ts); % 计算系数 A
% 仿真输入信号 (正弦和噪声信号的叠加)
t = 0:Ts:1-Ts; % 离散时间向量
input_signal = sin(2*pi*20*t) + sin(2*pi*100*t) + 0.5*randn(size(t)); % 混合信号: 5 Hz + 50 Hz + 噪声
% 初始化输出信号
output_signal = zeros(size(input_signal));
% 初始条件设置
output_signal(1) = input_signal(1); % 假设初始输出等于输入
% 递推公式计算一阶高通滤波器输出
for k = 2:length(input_signal)
output_signal(k) = A * output_signal(k-1) + A * (input_signal(k) - input_signal(k-1));
end
% FFT 分析
n = length(t); % 数据点数
f = (0:n-1)*(fs/n); % 频率向量
input_fft = abs(fft(input_signal)); % 输入信号的 FFT
output_fft = abs(fft(output_signal)); % 输出信号的 FFT
% 画图 - 输入输出信号的时域响应
figure;
subplot(2, 1, 1);
plot(t, input_signal);
title('输入信号 (时域)');
xlabel('时间 (s)');
ylabel('振幅');
subplot(2, 1, 2);
plot(t, output_signal);
title('滤波后信号 (时域)');
xlabel('时间 (s)');
ylabel('振幅');
% 画图 - 输入输出信号的频域响应
figure;
subplot(2, 1, 1);
plot(f(1:n/2), input_fft(1:n/2)); % 只显示正频率部分
title('输入信号 (频域)');
xlabel('频率 (Hz)');
ylabel('幅度');
subplot(2, 1, 2);
plot(f(1:n/2), output_fft(1:n/2)); % 只显示正频率部分
title('滤波后信号 (频域)');
xlabel('频率 (Hz)');
ylabel('幅度');
高通滤波器后信号的低频部分被显著抑制,只保留了高频成分。
通过 FFT 分析,可以看到高频成分(如 100 Hz)被保留,而低频成分(如 20 Hz)被滤除。
3、二阶低通滤波器
3.1 二阶RC低通滤波器的连续域数学模型
3.2 二阶RC低通滤波器的算法推导
3.3 matlab仿真
% 参数设置
fs = 1000; % 采样频率 1000 Hz
Ts = 1/fs; % 采样周期
w0 = 2 * pi * 50; % 截止频率 10 Hz
xi = 0.707; % 阻尼系数 (一般取0.707, 即1/sqrt(2))
% 计算滤波器系数
b0 = w0^2 * Ts^2;
a0 = w0^2 * Ts^2 + 4 * xi * w0 * Ts + 4;
a1 = 2 * w0^2 * Ts^2 - 8;
a2 = w0^2 * Ts^2 - 4 * xi * w0 * Ts + 4;
% 仿真输入信号 (正弦和噪声信号的叠加)
t = 0:Ts:1-Ts; % 离散时间向量
input_signal = sin(2*pi*20*t) + sin(2*pi*100*t) + 0.5*randn(size(t)); % 混合信号: 5 Hz + 50 Hz + 噪声
% 初始化输出信号
output_signal = zeros(size(input_signal));
% 初始条件设置
output_signal(1) = input_signal(1); % 假设初始输出等于输入
output_signal(2) = input_signal(2);
% 递推公式计算二阶低通滤波器输出
for k = 3:length(input_signal)
output_signal(k) = (b0/a0) * input_signal(k) + (2*b0/a0) * input_signal(k-1) + (b0/a0) * input_signal(k-2) ...
- (a1/a0) * output_signal(k-1) - (a2/a0) * output_signal(k-2);
end
% FFT 分析
n = length(t); % 数据点数
f = (0:n-1)*(fs/n); % 频率向量
input_fft = abs(fft(input_signal)); % 输入信号的 FFT
output_fft = abs(fft(output_signal)); % 输出信号的 FFT
% 画图 - 输入输出信号的时域响应
figure;
subplot(2, 1, 1);
plot(t, input_signal);
title('输入信号 (时域)');
xlabel('时间 (s)');
ylabel('振幅');
subplot(2, 1, 2);
plot(t, output_signal);
title('滤波后信号 (时域)');
xlabel('时间 (s)');
ylabel('振幅');
% 画图 - 输入输出信号的频域响应
figure;
subplot(2, 1, 1);
plot(f(1:n/2), input_fft(1:n/2)); % 只显示正频率部分
title('输入信号 (频域)');
xlabel('频率 (Hz)');
ylabel('幅度');
subplot(2, 1, 2);
plot(f(1:n/2), output_fft(1:n/2)); % 只显示正频率部分
title('滤波后信号 (频域)');
xlabel('频率 (Hz)');
ylabel('幅度');
滤波后信号的时域响应相比于输入信号的噪声平滑很多 , 通过 FFT 分析可以观察到 100 Hz 的高频成分被衰减,而 20 Hz 的低频成分保留。
4、二阶高通滤波器
4.1 二阶RC高通滤波器的连续域数学模型
4.2 二阶RC高通滤波器的算法推导
4.3 matlab仿真
% 参数设置
fs = 1000; % 采样频率 1000 Hz
Ts = 1/fs; % 采样周期
w0 = 2 * pi * 50; % 截止频率 10 Hz
xi = 0.707; % 阻尼系数 (一般取0.707, 即1/sqrt(2))
% 计算滤波器系数
b0 = 4;
a0 = w0^2 * Ts^2 + 4 * xi * w0 * Ts + 4;
a1 = 2 * w0^2 * Ts^2 - 8;
a2 = w0^2 * Ts^2 - 4 * xi * w0 * Ts + 4;
% 仿真输入信号 (正弦和噪声信号的叠加)
t = 0:Ts:1-Ts; % 离散时间向量
input_signal = sin(2*pi*20*t) + sin(2*pi*100*t) + 0.5*randn(size(t)); % 混合信号: 5 Hz + 50 Hz + 噪声
% 初始化输出信号
output_signal = zeros(size(input_signal));
% 初始条件设置
output_signal(1) = input_signal(1); % 假设初始输出等于输入
output_signal(2) = input_signal(2);
% 递推公式计算二阶高通滤波器输出
for k = 3:length(input_signal)
output_signal(k) = (b0/a0) * input_signal(k) - (2*b0/a0) * input_signal(k-1) + (b0/a0) * input_signal(k-2) ...
- (a1/a0) * output_signal(k-1) - (a2/a0) * output_signal(k-2);
end
% FFT 分析
n = length(t); % 数据点数
f = (0:n-1)*(fs/n); % 频率向量
input_fft = abs(fft(input_signal)); % 输入信号的 FFT
output_fft = abs(fft(output_signal)); % 输出信号的 FFT
% 画图 - 输入输出信号的时域响应
figure;
subplot(2, 1, 1);
plot(t, input_signal);
title('输入信号 (时域)');
xlabel('时间 (s)');
ylabel('振幅');
subplot(2, 1, 2);
plot(t, output_signal);
title('滤波后信号 (时域)');
xlabel('时间 (s)');
ylabel('振幅');
% 画图 - 输入输出信号的频域响应
figure;
subplot(2, 1, 1);
plot(f(1:n/2), input_fft(1:n/2)); % 只显示正频率部分
title('输入信号 (频域)');
xlabel('频率 (Hz)');
ylabel('幅度');
subplot(2, 1, 2);
plot(f(1:n/2), output_fft(1:n/2)); % 只显示正频率部分
title('滤波后信号 (频域)');
xlabel('频率 (Hz)');
ylabel('幅度');
滤波后信号的时域响应中,低频部分被显著衰减 ,通过 FFT 分析,可以观察到高频部分的保留以及低频部分的衰减 。
5、陷波滤波
陷波滤波器其实就是带阻滤波器。
这里只提供matlab算法验证,详细公式推导前看《陷波滤波器的数学模型推导及算法实现》
% 参数设置
fs = 1000; % 采样频率 1000 Hz
Ts = 1/fs; % 采样周期
w_n = 2 * pi * 50; % 陷波频率 50 Hz
B = 2 * pi * 20; % 带宽 (单位: rad)
depth = 0.05; % 深度 (在范围 -sqrt(2)/2 到 sqrt(2)/2 之间) 越小效果越好
% 计算 k_1 和 k_2
k1 = sqrt((1 - sqrt(1 + (B^2 / w_n^2))) / (4 * (depth)^2 - 2));
k2 = k1 * depth;
% 计算滤波器系数
a_1 = w_n^2 * Ts^2 + 4 * Ts * k1 * w_n + 4;
a_2 = 2 * w_n^2 * Ts^2 - 8;
a_3 = w_n^2 * Ts^2 - 4 * Ts * k1 * w_n + 4;
b_1 = w_n^2 * Ts^2 + 4 * Ts * k2 * w_n + 4;
b_2 = 2 * w_n^2 * Ts^2 - 8;
b_3 = w_n^2 * Ts^2 - 4 * Ts * k2 * w_n + 4;
% 生成测试信号 (正弦和噪声信号的叠加)
t = 0:Ts:1-Ts; % 离散时间向量
input_signal = sin(2*pi*20*t) +sin(2*pi*50*t)+ sin(2*pi*100*t) + 0.5*randn(size(t)); % 混合信号: 20 Hz + 100 Hz + 噪声
% 初始化输出信号
output_signal = zeros(size(input_signal));
% 初始条件设置
output_signal(1) = input_signal(1);
output_signal(2) = input_signal(2);
% 递推公式计算陷波滤波器输出
for n = 3:length(input_signal)
output_signal(n) = (b_1/a_1) * input_signal(n) + (b_2/a_1) * input_signal(n-1) + (b_3/a_1) * input_signal(n-2) ...
- (a_2/a_1) * output_signal(n-1) - (a_3/a_1) * output_signal(n-2);
end
% FFT 分析
n = length(t); % 数据点数
f = (0:n-1)*(fs/n); % 频率向量
input_fft = abs(fft(input_signal)); % 输入信号的 FFT
output_fft = abs(fft(output_signal)); % 输出信号的 FFT
% 画图 - 输入输出信号的时域响应
figure;
subplot(2, 1, 1);
plot(t, input_signal);
title('输入信号 (时域)');
xlabel('时间 (s)');
ylabel('振幅');
subplot(2, 1, 2);
plot(t, output_signal);
title('滤波后信号 (时域)');
xlabel('时间 (s)');
ylabel('振幅');
% 画图 - 输入输出信号的频域响应
figure;
subplot(2, 1, 1);
plot(f(1:n/2), input_fft(1:n/2)); % 只显示正频率部分
title('输入信号 (频域)');
xlabel('频率 (Hz)');
ylabel('幅度');
subplot(2, 1, 2);
plot(f(1:n/2), output_fft(1:n/2)); % 只显示正频率部分
title('滤波后信号 (频域)');
xlabel('频率 (Hz)');
ylabel('幅度');
通过 FFT 分析,可以观察到50Hz部分的衰减,而20Hz和100Hz得到了保留。
6、带通滤波器
6.1 带通滤波器的连续域数学模型
6.2 带通滤波器的算法推导
6.3 matlab仿真
% 参数设置
fs = 1000; % 采样频率 1000 Hz
Ts = 1/fs; % 采样周期
f_0 = 70; % 中心频率 50 Hz
BW = 2 * pi * 20; % 带宽 20 Hz
k = 0.5; % 增益常数
% 计算角频率和品质因数
omega_0 = 2 * pi * f_0;
Q = omega_0 / BW;
% 计算数字滤波器系数
b_0 = 2 * Ts * k * BW;
a_2 = 4 - 2 * BW * Ts + omega_0^2 * Ts^2;
a_1 = 2 * Ts^2 * omega_0^2 - 8;
a_0 = 4 + 2 * BW * Ts + Ts^2 * omega_0^2;
% 生成测试信号 (正弦和噪声信号的叠加)
t = 0:Ts:1-Ts; % 离散时间向量
input_signal = sin(2*pi*20*t) + sin(2*pi*70*t) + sin(2*pi*120*t) + 0.5*randn(size(t)); % 混合信号: 20 Hz + 50 Hz + 100 Hz + 噪声
% 初始化输出信号
output_signal = zeros(size(input_signal));
% 初始条件设置
output_signal(1) = input_signal(1);
output_signal(2) = input_signal(2);
% 递推公式计算带通滤波器输出
for n = 3:length(input_signal)
output_signal(n) = (b_0/a_0) * input_signal(n) + (b_0/a_0) * input_signal(n-2) ...
- (a_1/a_0) * output_signal(n-1) - (a_2/a_0) * output_signal(n-2);
end
% FFT 分析
n = length(t); % 数据点数
f = (0:n-1)*(fs/n); % 频率向量
input_fft = abs(fft(input_signal)); % 输入信号的 FFT
output_fft = abs(fft(output_signal)); % 输出信号的 FFT
% 画图 - 输入输出信号的时域响应
figure;
subplot(2, 1, 1);
plot(t, input_signal);
title('输入信号 (时域)');
xlabel('时间 (s)');
ylabel('振幅');
subplot(2, 1, 2);
plot(t, output_signal);
title('滤波后信号 (时域)');
xlabel('时间 (s)');
ylabel('振幅');
% 画图 - 输入输出信号的频域响应
figure;
subplot(2, 1, 1);
plot(f(1:floor(n/2)), input_fft(1:floor(n/2))); % 只显示正频率部分
title('输入信号 (频域)');
xlabel('频率 (Hz)');
ylabel('幅度');
subplot(2, 1, 2);
plot(f(1:floor(n/2)), output_fft(1:floor(n/2))); % 只显示正频率部分
title('滤波后信号 (频域)');
xlabel('频率 (Hz)');
ylabel('幅度');
给定信号是20Hz、70Hz、120Hz以及噪声信号组成,带通滤波器中心频率为70Hz,经过带通滤波器后的信号通过FFT分析对比,20Hz和120Hz的信号进行了衰减,保留了70Hz信号。
参考
【1】基于STM32的ADC采样及各式滤波实现(HAL库,含VOFA+教程):
https://blog.csdn.net/black_sneak/article/details/129629485
【2】【学习笔记】matlab进行数字信号处理(三)数字滤波技术:https://blog.csdn.net/weixin_42853410/article/details/114407188?share_token=2D1ED73C-F87F-4B9D-BCEF-9AFEE203C6C9&tt_from=copy_link&utm_source=copy_link&utm_medium=toutiao_ios&utm_campaign=client_share 【学习笔记】matlab进行数字信号处理(三)数字滤波技术_信号滤波技术 - 今日头条
【3】RC低通滤波器截止频率公式推导:
https://blog.csdn.net/zwc475793240/article/details/122432824
【4】一阶RC高通滤波器详解(仿真+matlab+C语言实现):
https://bbs.huaweicloud.com/blogs/314645
【5】滤波器——二阶低通滤波器:
https://blog.csdn.net/qq_53043199/article/details/136071435#:~:text=%E6%9C%AC%E6%96%87%E8%AF%A6%E7%BB%86%E4%BB%8B%E7%BB%8D%E4%BA%86%E4%BA%8C
【6】数字二阶低通滤波器公式推导及代码实现
https://blog.csdn.net/qq_26988431/article/details/100779047#:~:text=%E4%BA%8C%E9%98%B6%E4%BD%8E%E9%80%9A%E6%BB%A4%E6%B3%A2%E5%99%A8%E5%8F%AF
【7】陷波滤波器(Notch Filter)的离散化设计:
https://blog.csdn.net/u013581448/article/details/116743786?utm_medium=distribute.pc_relevant.none-task-blog-2defaultbaidujs_baidulandingword~default-4-116743786-blog-105431862.235v38pc_relevant_anti_t3&spm=1001.2101.3001.4242.3&utm_relevant_index=7
【8】A Better Way to Think About a Notch Filter | Control Systems in Practice
https://www.youtube.com/watch?v=tpAA5eUb6eo