二阶巴特沃斯低通滤波器

目录引言模拟二阶巴特沃斯低通滤波器双线性变换的基本公式为什么需要预畸变预畸变的直观理解将模拟传递函数代入双线性变换归一化成标准数字滤波器形式更常用的简化形式最终的差分方程巴特沃斯滤波器完整设计流程Matlab分析验证幅频响应图上图相频响应图下图延时补偿一阶段预测补偿MCU代码实现测试结果引言巴特沃斯滤波器是一种模拟滤波器它在频率响应方面具有特殊的属性。它被设计为具有均匀的幅度响应即在通带内它对所有频率的增益是相等的而在阻带内它对所有频率的增益都是零。为了实现这种频率响应巴特沃斯滤波器在设计时采用了一些特殊的方法。具体而言它是由一组极点和零点组成的这些极点和零点的位置是根据所需的通带和阻带响应来计算的。在巴特沃斯滤波器中极点和零点都是共轭对这意味着它们的实部相等虚部互为相反数。巴特沃斯滤波器被广泛用于模拟信号处理中例如在音频处理中进行低通滤波。但是现代数字信号处理技术已经使得数字滤波器成为实际上的选择因为数字滤波器可以实现更复杂的频率响应而且具有更好的可编程性和可重复性。模拟二阶巴特沃斯低通滤波器给定模拟域传递函数其中模拟截止角频率单位是。二阶巴特沃斯低通的特点是也就是在截止频率处衰减。双线性变换的基本公式双线性变换把模拟域 s平面映射到数字域 z平面其中所以也可以写为为了后面推导方便定义于是为什么需要预畸变双线性变换有一个很重要的问题它会把模拟频率轴 Ω非线性压缩到数字频率轴 ω。模拟角频率和数字角频率之间的关系是也就是其中Ω模拟角频率单位 rad/s数字角频率单位 rad/sample如果数字滤波器希望截止频率为代入频率映射关系中我们可以得到因此这就是预畸变频率。预畸变的直观理解你想让数字滤波器在处刚好是。但是双线性变换会扭曲频率轴。所以不能简单令而应该先算一个“修正后的模拟截止频率”然后把这个 Ωc放进模拟传递函数这样经过双线性变换后数字滤波器的截止频率才会准确落在。当时所以也就是说低频情况下预畸变影响不大但当截止频率接近采样频率的较高比例时预畸变就非常重要。将模拟传递函数代入双线性变换我们定义双线性变换就可以写成模拟滤波器为代入 s为了整理成多项式形式分子分母同乘得到展开分子展开也就是展开分母第一项第二项第三项定义变量:对这些变量进行合并:最终我们就可以得到归一化成标准数字滤波器形式数字滤波器一般形式为所以我们给式子的分子分母同时除以D0,我们便可以得到对应的一个系数:其中更常用的简化形式我们定义因为而所以将所有系数除以可得于是二阶巴特沃斯低通数字滤波器系数为其中最终的差分方程由我们可以得到我们将其转换到时间域所以巴特沃斯滤波器完整设计流程首先需要确定第一步计算归一化预畸变参数第二步计算公共分母第三步计算系数第四步代入差分方程Matlab分析验证对上述推导进行验证 选择采样频率为20khz 截至频率为200hz 滤波器阶数为2阶%% 二阶巴特沃斯低通滤波器设计与验证 % 采样频率 fs 20 kHz % 截止频率 fc 200 Hz % 方法模拟巴特沃斯原型 预畸变 双线性变换 clear; clc; close all; %% 1. 基本参数 fs 20000; % 采样频率Hz fc 200; % 截止频率Hz Ts 1 / fs; % 采样周期s N 2; % 二阶巴特沃斯滤波器 %% 2. 预畸变处理 % 双线性变换会造成频率轴压缩因此需要预畸变 % 目标保证数字滤波器在 fc 200 Hz 处刚好为 -3 dB Omega_c 2 * fs * tan(pi * fc / fs); % 预畸变后的模拟截止角频率rad/s fprintf(预畸变后的模拟截止角频率 Omega_c %.6f rad/s\n, Omega_c); fprintf(未预畸变时 2*pi*fc %.6f rad/s\n, 2*pi*fc); %% 3. 生成模拟域二阶巴特沃斯低通滤波器 % butter(N, Wn, s) 中 Wn 的单位是 rad/s [bs, as] butter(N, Omega_c, s); disp(模拟域滤波器分子系数 bs ); disp(bs); disp(模拟域滤波器分母系数 as ); disp(as); %% 4. 使用双线性变换离散化 % 注意这里已经手动完成了预畸变因此 bilinear 中不要再次传入预畸变频率 [bz, az] bilinear(bs, as, fs); disp(数字滤波器分子系数 bz ); disp(bz); disp(数字滤波器分母系数 az ); disp(az); %% 5. 打印差分方程系数 b0 bz(1); b1 bz(2); b2 bz(3); a1 az(2); a2 az(3); fprintf(\n数字滤波器差分方程为\n); fprintf(y[n] %.12f*x[n] %.12f*x[n-1] %.12f*x[n-2] , b0, b1, b2); fprintf(- (%.12f)*y[n-1] - (%.12f)*y[n-2]\n, a1, a2); fprintf(\n标准系数\n); fprintf(b0 %.12f\n, b0); fprintf(b1 %.12f\n, b1); fprintf(b2 %.12f\n, b2); fprintf(a1 %.12f\n, a1); fprintf(a2 %.12f\n, a2); %% 6. 频域响应分析 nfft 8192; [H, f] freqz(bz, az, nfft, fs); mag_dB 20 * log10(abs(H)); phase_deg unwrap(angle(H)) * 180/pi; % 找到最接近截止频率 fc 的频点 [~, idx_fc] min(abs(f - fc)); mag_fc mag_dB(idx_fc); phase_fc phase_deg(idx_fc); fprintf(\n截止频率 %.2f Hz 处\n, fc); fprintf(幅值 %.6f dB\n, mag_fc); fprintf(相位 %.6f deg\n, phase_fc); %% 7. 绘制幅频响应 figure(Name, Magnitude Response); plot(f, mag_dB, LineWidth, 1.5); grid on; xlabel(Frequency / Hz); ylabel(Magnitude / dB); title(Second-Order Butterworth Low-Pass Filter Magnitude Response); xlim([0 2000]); ylim([-60 5]); hold on; plot(fc, mag_fc, ro, MarkerSize, 8, LineWidth, 1.5); xline(fc, --r, Cutoff Frequency); yline(-3, --k, -3 dB); text(fc 40, mag_fc, ... sprintf(fc %.0f Hz\\nMagnitude %.3f dB, fc, mag_fc), ... FontSize, 10); hold off; %% 8. 绘制相频响应 figure(Name, Phase Response); plot(f, phase_deg, LineWidth, 1.5); grid on; xlabel(Frequency / Hz); ylabel(Phase / degree); title(Second-Order Butterworth Low-Pass Filter Phase Response); xlim([0 2000]); hold on; plot(fc, phase_fc, ro, MarkerSize, 8, LineWidth, 1.5); xline(fc, --r, Cutoff Frequency); text(fc 40, phase_fc, ... sprintf(fc %.0f Hz\\nPhase %.3f deg, fc, phase_fc), ... FontSize, 10); hold off; %% 9. 极零点分析 % 使用 zplane 绘制零极点图 figure(Name, Pole-Zero Plot by zplane); zplane(bz, az); grid on; title(Pole-Zero Plot in z-plane); %% 10. 使用 pzmap 绘制零极点图 % 需要 Control System Toolbox sysz tf(bz, az, Ts, Variable, z^-1); figure(Name, Pole-Zero Plot by pzmap); pzmap(sysz); grid on; title(Pole-Zero Map of Digital Butterworth Filter); %% 11. 稳定性判断 poles_z roots(az); zeros_z roots(bz); disp(数字滤波器零点); disp(zeros_z); disp(数字滤波器极点); disp(poles_z); pole_abs abs(poles_z); disp(极点模值); disp(pole_abs); if all(pole_abs 1) disp(稳定性判断所有极点均位于单位圆内该数字滤波器稳定。); else disp(稳定性判断存在极点不在单位圆内该数字滤波器不稳定。); end %% 12. 与理论 -3 dB 对比 fprintf(\n理论上二阶巴特沃斯滤波器在截止频率处应为 -3.0103 dB。\n); fprintf(当前仿真结果为 %.6f dB。\n, mag_fc); fprintf(误差为 %.6f dB。\n, mag_fc - (-3.0103));预畸变后的模拟截止角频率 Omega_c 1257.050642 rad/s 未预畸变时 2*pi*fc 1256.637061 rad/s 模拟域滤波器分子系数 bs 1.0e06 * 0 0 1.5802 模拟域滤波器分母系数 as 1.0e06 * 0.0000 0.0018 1.5802 数字滤波器分子系数 bz 0.0009 0.0019 0.0009 数字滤波器分母系数 az 1.0000 -1.9112 0.9150 数字滤波器差分方程为 y[n] 0.000944691844*x[n] 0.001889383688*x[n-1] 0.000944691844*x[n-2] - (-1.911197067426)*y[n-1] - (0.914975834801)*y[n-2] 标准系数 b0 0.000944691844 b1 0.001889383688 b2 0.000944691844 a1 -1.911197067426 a2 0.914975834801 截止频率 200.00 Hz 处 幅值 -3.018792 dB 相位 -90.079143 deg 数字滤波器零点 -1.0000 0.0000i -1.0000 - 0.0000i 数字滤波器极点 0.9556 0.0425i 0.9556 - 0.0425i 极点模值 0.9565 0.9565 稳定性判断所有极点均位于单位圆内该数字滤波器稳定。 理论上二阶巴特沃斯滤波器在截止频率处应为 -3.0103 dB。 当前仿真结果为 -3.018792 dB。 误差为 -0.008492 dB。幅频响应图上图通带特性0~200Hz截止频率 fc200Hz为通带低频信号增益接近 0dB几乎无衰减信号完整通过截止点标准特性200Hz 处幅值约 - 3.01dB完美符合巴特沃斯滤波器 - 3dB 截止定义阻带衰减频率超过 200Hz 后幅值快速向下衰减频率越高衰减越剧烈高频分量被大幅抑制二阶特征曲线斜率对应二阶滤波高频处每倍频程衰减约 - 40dB / 十倍频滤波抑制能力优于一阶低通。相频响应图下图相位滞后规律输入信号经过滤波器后会产生相位滞后频率越高滞后越严重相位数值越负截止频率相位200Hz 截止点相位为 -90.079°极限相位当频率远大于截止频率时相位趋近 - 180°非线性相位相位曲线为弯曲曲线不是直线代表群延迟随频率变化—— 不同频率信号延时不一样会造成波形相位失真在电机 FOC、同步采样场景需要重点关注。延时补偿二阶巴特沃斯低通在低频段的延时近似为假如截止频率是那么:如果采样频率是采样周期就是那么延时大概是也就是说你这个二阶低通大约会带来20 多个采样周期的相位滞后。一阶段预测补偿在工程上最简单、最适合 MCU 的方式是也可以写成其中这个公式的含义是低通滤波后的信号已经滞后了所以用它当前的变化趋势向前预测一点。MCU代码实现/***********h文件*************/ typedef struct { float Raw_Sign; //输入原始信号 float Filter_Sign; //滤波输出信号 float Comp_Sign; //延时补偿信号 /*滤波器系数*/ float b0; float b1; float b2; float a1; float a2; /*中间保存变量*/ float z1; float z2; float y_last; //上一次滤波输出 float Delay_comp_k; //延时补偿系数 float fc; //截至频率 float fs; //采样频率 }Butterwoth_LPF2; //二阶巴特沃斯滤波器 /***********c文件*************/ //巴特沃斯滤波器系数计算 void Foc_Butterworth_LPF2_CallCoffs(Butterwoth_LPF2 *BLp) { float Theta; //角度中间变量 float Theta2; float D; float Td; float Ts; //计算公用系数 Theta tanf(PI * BLp-fc / BLp-fs); Theta2 Theta * Theta; D 1.0f sqrt2 * Theta Theta2; BLp-b0 Theta2 / D; BLp-b1 2.0f * Theta2 / D; BLp-b2 Theta2 / D; BLp-a1 2.0f * (Theta2 - 1.0f) / D; BLp-a2 (1.0f - sqrt2 * Theta Theta2)/D; //补偿系数计算 Ts 1.0f / BLp-fs; Td sqrt2 / (2.0f * PI *BLp-fc); BLp-Delay_comp_k Td/Ts; } //巴特沃斯滤波器参数全局初始化 void Foc_Butterworth_LPF2_Init(Butterwoth_LPF2 *BLp) { BLp-Raw_Sign 0.0f; BLp-Filter_Sign 0.0f; BLp-Comp_Sign 0.0f; BLp-b0 0.0f; BLp-b1 0.0f; BLp-b2 0.0f; BLp-a1 0.0f; BLp-a2 0.0f; BLp-z1 0.0f; BLp-z2 0.0f; BLp-y_last 0.0f; Foc_Butterworth_LPF2_CallCoffs(BLp); } //滤波器数据更新 void Foc_Butterworth_LPF2_Update(Butterwoth_LPF2 *BLp) { float y; float y_comp; y BLp-b0 * BLp-Raw_Sign BLp-z1; BLp-z1 BLp-b1 * BLp-Raw_Sign - BLp-a1 * y BLp-z2; BLp-z2 BLp-b2 * BLp-Raw_Sign - BLp-a2 *y; BLp-Filter_Sign y; /*一阶延时补偿*/ y_comp y BLp-Delay_comp_k * (y - BLp-y_last); BLp-y_last y; BLp-Comp_Sign y_comp; }测试结果本次测试使用mcu为瑞萨R7A6系列主控 生成原始信号为幅值为1频率为50hz的正弦波生成的高频噪声为幅值为0.3频率为2000hz的正弦波。滤波器截至频率为50hz采样执行频率为20000hz满足奈奎斯特采样定理。50hz 原始信号叠加2000hz的信号滤波后的信号和原始信号存在1ms相位延迟经过时间补偿后的原始信号和滤波信号相位延时已经消除觉得有用的话帮忙点个赞啦谢谢大家了