自学内容网 自学内容网
自学内容网 > 正文

二阶滤波算法总结（对RC滤波算法整理的部分修正和完善）

🕗 发布于 2024-09-25 20:24 算法  滤波  信号处理  控制  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+wc​Ts​wc​Ts​​
上式就是已知截至频率 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+wc​Ts​1​
上式就是已知截至频率 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


原文地址：https://blog.csdn.net/qq_28149763/article/details/142499020

免责声明：本站文章内容转载自网络资源，如本站内容侵犯了原著者的合法权益，可联系本站删除。更多内容请关注自学内容网（zxcms.com）！

相关文章

自学内容网
Copyright © 2015-2024