简介《课题二 基于MATLAB平台的心电信号分析系统设计及仿真》是一份面向电子信息、生物医学工程等专业课程设计与毕业设计的PDF技术文档适合正在学习数字信号处理、MATLAB编程及Simulink仿真的学生与工程人员参考。文档围绕MIT-BIH心电数据库展开完整设计流程从txt格式心电数据的读取与实际波形还原到线性插值处理使非均匀采样数据等距化再依据心电信号频域特点设计低通、高通及带通滤波器并结合巴特沃斯或切比雪夫滤波器分析幅频特性、零极点分布与阶跃响应。文中还给出50Hz工频陷波器的设计思路以及直接编程静态仿真与Simulink动态建模两条实现路径最后对处理前后的信号做频谱对比分析并整理课程设计报告撰写要求与参考文献。资源包内仅含1个PDF文件大小约1019KB便于在电脑或移动端直接查阅。目前已有357人学习适合希望系统掌握心电信号预处理、滤波器设计与仿真验证方法的读者按步骤实践。1. 心电信号不是丢进滤波器就能出结果MIT-BIH 的 txt 文件长得挺朴素第一列采样时刻第二列 MLII 导联幅值第三列 V5 导联幅值Tab 分隔。坑不在读取而在时间轴——它不是等间距的采集丢点、标注跳变都会让 Δt 在 0.8ms 到 3ms 之间漂移。而 butter、cheby1、ellip 这些函数默认输入是按 Fs 均匀采样的序列你拿非均匀数据直接滤通带边缘会被抹平ST 段那几十微伏的细节直接淹没在数值误差里。所以这个课题真正的执行顺序是读数据并跳过说明行 → 一次线性插值成 1ms 等间隔 → 设计模拟原型滤波器 → 双线性变换离散化 → 50Hz 陷波 → Simulink 里跑动态仿真 → 用 FFT 前后对比频谱反推参数。适合做信号与系统课程设计、DSP 大作业以及刚上手 MATLAB 信号处理工具箱的人。2. MIT-BIH txt 读取与 1ms 等间隔线性插值2.1 三列数据读取与说明行跳过原始文件前两行是文字说明不是数值。要求是严禁手动删除必须用程序忽略。老教材常用textread但它在新版本里已经不推荐我一般直接用textscan配fgetl。fid fopen(100.txt, r); for k 1:2 fgetl(fid); % 只读走前两行说明文字不解析 end C textscan(fid, %f %f %f); % 三列时间、MLII、V5 fclose(fid); t C{1}; % 采样时刻单位秒 x C{2}; % MLII 导联幅值 % C{3} 是 V5同一心拍的另一导联分析只留一路即可textscan返回的是 cell 数组所以要用C{1}取列。如果你把%f %f %f写成%f%f%f遇到 Tab 分隔同样能解析但列数对不上时会静默错位建议保留空格。读完后先做一次体检any(isnan(x))如果返回 1基本就是说明行没跳干净或者文件里混了空行。常见做法是再补一句x x(~isnan(x)); t t(~isnan(t));兜底。2.2 一次线性插值的两种实现插值公式不复杂相邻两点时间差 Δt需要的插值点数 N Δt/0.001幅值步进 ΔA/N然后从 tj-1 ti-1、Aj-1 Ai-1 逐点推。手写循环更能体现原理也方便答辩时讲清楚。Ts 0.001; % 目标采样间隔 1ms tq t(1); xq x(1); % 预置首点 for i 2:length(t) dt t(i) - t(i-1); N round(dt / Ts); % 两点之间需要插几个点 if N 1 % 已经够密不插直接接上 tq(end1,1) t(i); xq(end1,1) x(i); continue; end for j 1:N-1 tq(end1,1) t(i-1) j*Ts; xq(end1,1) x(i-1) (x(i)-x(i-1))*j/N; end tq(end1,1) t(i); xq(end1,1) x(i); endround而不是floor是关键用floor会让累计时间轴越走越慢10 秒数据末尾能偏出几十毫秒。上面这种end1追加写法在大数据量下效率很差正式跑之前把tq、xq预分配成zeros(sum(round(diff(t)/Ts))1, 1)会快一个量级。工程里更省事的是直接调interp1tq t(1):Ts:t(end); xq interp1(t, x, tq, linear); % 默认就是 linear写出来更清楚两者结果在数值上一致差别只是interp1对t中重复时刻点会报警告而手写循环会算出 N0 然后跳过。如果数据里确实有时间戳重复先[t, idx] unique(t); x x(idx);去重再插。2.3 插值结果的校验插完不要急着画图交差。先验三件事max(diff(tq))应该等于 0.001浮点误差范围内如果出现 0.002说明有段原始数据 Δt 小于 0.5ms 被round吞掉了。length(tq)约等于t(end)/0.0013600 点原始数据大约插到 10000 点上下。corr(x, interp1(tq, xq, t))接近 1说明插值没有引入明显畸变。fprintf(插值后点数 %d最大间隔 %.6f s\n, numel(tq), max(diff(tq))); plot(t, x, k., tq, xq, r-); % 原始点叠加插值曲线 xlabel(时间 / s); ylabel(幅值 / mV); legend(原始,插值后);把插值后的[tq, xq]另存为ecg_interp.txt用dlmwrite(ecg_interp.txt, [tq xq], delimiter, \t, precision, %.6f)后面 Simulink 和滤波程序都从这里读避免每次重复插值。3. 巴特沃斯与切比雪夫原型及双线性变换3.1 从心电频谱倒推通带阻带指标正常人心电能量集中在 0.05Hz 到 100HzQRS 波主频在 1025HzP 波和 T 波更低。幅度上成人约 1mV胎儿只有 10μV 量级。这意味着两件事低频端要压住呼吸引起的 0.150.5Hz 基线漂移高频端要挡掉肌电和工频。常见指标组合是低通 Fpass100Hz、Fstop102Hz高通 Fpass1Hz、Fstop0.5Hz阻带衰减 40dB 以上。MIT-BIH 的采样率是 360Hz奈奎斯特 180Hz100Hz 通带留了足够过渡带。滤波器类型通带/阻带阶数作用高通巴特沃斯 IIRFpass1Hz, Fstop0.5Hz4去基线漂移低通等波纹 FIRFpass100Hz, Fstop102Hz最小阶去高频肌电陷波巴特沃斯带阻4853Hz4去 50Hz 工频选巴特沃斯是因为通带最平坦心电的 ST 段对通带波动很敏感1dB 的起伏就能让抬高/压低判断出错。切比雪夫 I 型过渡带更陡同样阶数下能压得更狠代价是通带等波纹用在陷波这种窄带场景更划算。3.2 模拟原型与 freqs 幅频特性先在模拟域把原型搭出来画freqs曲线验证指标这一步是答辩里最容易加分的地方。Rp 1; Rs 40; % 通带 1dB阻带 40dB Wp 2*pi*100; Ws 2*pi*130; % 模拟角频率 rad/s [N, Wc] buttord(Wp, Ws, Rp, Rs, s); [bs, as] butter(N, Wc, s); % 模拟低通原型系数 w logspace(0, 3, 800); % 1 rad/s 到 1000 rad/s H freqs(bs, as, w); semilogx(w/(2*pi), 20*log10(abs(H))); grid on; xlabel(频率 / Hz); ylabel(幅度 / dB);buttord里带s才返回模拟域阶数和截止频率不带就是数字域。Wc通常落在Wp和Ws之间别拿它当 Fpass 用。切比雪夫版本只需换成[N, Wc] cheb1ord(Wp, Ws, Rp, Rs, s)和cheby1(N, Rp, Wc, s)注意cheby1要显式传通带波纹 Rp。freqs给的是模拟频率响应看的是原型形状对不对。零极点分布用[z, p, k] buttap(N)拿归一化原型再把p乘Wc去归一化。所有极点都应该落在 s 左半平面只要有一个跑到右半平面这个滤波器就是发散的不稳定系统必须回头查阶数或参数。3.3 双线性变换的预畸变模拟原型到数字滤波器有两条路。脉冲响应不变法时域逼近好但频谱会混叠只适合低通带通双线性变换把整个 s 平面压进单位圆没有混叠代价是频率轴非线性压缩Ω (2/T)·tan(ω/2)。所以数字截止频率必须先做预畸变再传给模拟原型。Fs 360; % MIT-BIH 采样率 fp 100; fsb 130; % 数字域 Hz Wp 2*Fs*tan(pi*fp/Fs); % 预畸变到模拟域 Ws 2*Fs*tan(pi*fsb/Fs); [N, Wc] buttord(Wp, Ws, Rp, Rs, s); [bs, as] butter(N, Wc, s); [bz, az] bilinear(bs, as, Fs); % 双线性离散化不预畸变直接bilinear的后果很隐蔽截止频率会往低频偏100Hz 的设计点实际可能落在 92Hz 附近通带边缘的心电高频成分被削掉QRS 幅度看起来变矮了。检查方法是freqz(bz, az, 1024, Fs)画数字幅频看 -3dB 点是不是还在 100Hz 附近。等效的写法是直接在数字域用buttord(Wp/(Fs/2), Ws/(Fs/2), Rp, Rs)再butterMATLAB 内部帮你做了预畸变结果和上面一致。两条路都可以走但报告里写清楚用了哪一种别混着放参数。4. Simulink 级联二阶节建模与 50Hz 陷波4.1 六阶低通拆成三个二阶节六阶 Butterworth 低通直接展开成一个六阶传递函数系数动态范围大定点实现或长时间仿真容易数值不稳定。标准做法是拆成三个二阶节级联A 0.09036; B [1.2686, 1.0106, 0.9044]; % 各节分母一次项系数 C [0.7051, 0.3583, 0.2155]; % 各节分母常数项 % 每节形式 Hk(z) A*(1 2z^-1 z^-2) / (1 - Bk*z^-1 Ck*z^-2)三个节串起来Numerator填[A 2*A A]Denominator填[1 -B(k) C(k)]。为什么分子统一是A(12z⁻¹z⁻²)因为 A 是总增益按节数均分后的结果(1z⁻¹)²保证每个节在 z-1 处有零点级联后在奈奎斯特频率处形成六重零点这正是 Butterworth 低通的典型特征。滤波器级数的顺序也有讲究。把 Q 值最高、极点最靠近单位圆的那一节放在最后能让前级先衰减掉一部分能量降低末级的动态范围压力。用zp2sos自动排序比自己拍脑袋排更靠谱[z, p, k] butter(N, Wn); sos zp2sos(z, p, k, down, up); % 按极点靠近单位圆程度排序4.2 From Workspace 与 FDATool 导入Simulink 里喂数据有两条路。一条是From Workspace模块变量名填simin结构体格式simin.time tq; simin.signals.values xq; simin.signals.dimensions 1;老版本要求是这种结构体新版本也接受timeseries(tq, xq)。如果报维度不匹配八成是values列向量写成了行向量加个(:)拉直。数据源记得设成0:Ts:t(end)的有限时长默认的inf会让仿真跑不完。另一条是Digital Filter Design模块先把 FDATool 里设计好的1.fda、2.fda、3.fda保存好在模块 Main 页的 Filter 栏选 Imported加载对应文件。FDATool 的参数对应关系是低通选 Lowpass FIR Equiripple Minimum Order填 Fpass100、Fstop102高通选 Highpass IIR Butterworth Minimum Order填 Fpass1、Fstop0.5带阻选 Bandstop IIR Butterworth填 Fstop149、Fpass148、Fpass253、Fstop251。注意FDATool 里 Minimum Order 算出来的阶数和你在命令行用 buttord 算出来的可能差 1因为 FDATool 用的是它自己的过渡带定义。报告里以 FDATool 导出的系数为准命令行版本只用来交叉验证。4.3 50Hz 陷波器的窄带挑战工频干扰是以 50Hz 为中心、带宽小于 1Hz 的窄带噪声。窄带带阻滤波器有个天然矛盾过渡带越窄阶数越高极点越贴近单位圆对系数量化越敏感。用 FDATool 设计 4853Hz 带阻但实际仿真发现 50Hz 只压下去 20dB 左右说明阶数不够。% 直接算 IIR 陷波比 fdatool 更可控 wo 50/(Fs/2); % 归一化中心频率 bw wo/35; % 3dB 带宽约 1.4Hz [b_n, a_n] iirnotch(wo, bw); y filtfilt(b_n, a_n, xq); % 零相位滤波避免 T 波被拖歪iirnotch的第二个参数是归一化带宽wo/35对应大约 1.4Hz改分母就能调窄。带宽小于 1Hz 时幅度响应会出现很深的凹口但相位响应在凹口附近剧烈变化用filter会造成 T 波形态失真所以这里必须用filtfilt双向滤波。代价是输出长度为原始两倍左右的延迟被抵消边界处会有3*(N-1)点的瞬态前后各截掉一段再算频谱。陷波前先确认 50Hz 干扰是不是真的存在。把原始信号做一次 FFT看 50Hz 处有没有明显尖峰。如果尖峰比周围噪声高不到 10dB先别设计陷波器多半是数据本身已经做过工频抑制硬加陷波反而会把 50Hz 附近的心电成分一起削掉。5. 频谱前后对比与参数回调做完滤波最终判断标准只有一个频谱对比。我一般把原始和滤波后的信号各做一次单边幅度谱用对数坐标叠在一起看。function [f, P] spec1(x, Fs) x x(:) - mean(x); % 去直流 N 2^nextpow2(length(x)); % 补到 2 的幂加快 FFT w hann(length(x)); % 加窗抑制泄漏 X fft(x.*w, N); P abs(X)/sum(w); % 幅度归一化 P P(1:N/21); P(2:end-1) 2*P(2:end-1); f (0:N/2)*(Fs/N); end [f1, P1] spec1(xq, 360); [f2, P2] spec1(y, 360); semilogy(f1, P1, k, f2, P2, r); xlabel(频率 / Hz); ylabel(幅度); legend(滤波前,滤波后);看三个位置就知道参数对不对。0.3Hz 附近滤波后应该下降 20dB 以上基线漂移被压住。50Hz凹口深度如果不到 30dB把iirnotch的带宽再收窄一档。100Hz 以上应该快速滚降如果还在缓降说明低通阶数不够。如果频谱看着没什么变化先别急着换滤波器类型按这个顺序回查时间轴单位。MIT-BIH 的 txt 里第一列单位是秒还是毫秒读错了 Fs 就是 0.36Hz频谱整个挤在左端。filtfilt的边界瞬态。开头和结尾各 200 点常常是数值爆炸区算频谱前先y y(200:end-200);。滤波器系数是否用了归一化频率。butter(N, 100)里的 100 是归一化到奈奎斯特的值必须是小于 1 的数写 100 会直接报错或返回全通。零极点位置。zplane(b, a)看极点有没有跑到单位圆外或者贴到 0.999 这种边缘位置后者说明阶数虚高实际已经临界稳定。最后给一个验证闭环的小技巧把滤波前后的信号各自做一次 QRS 检测哪怕只是最简单的阈值法比较 R 峰位置偏移。如果滤波后 R 峰整体平移超过 5ms说明滤波器群延迟没有被补偿重新用filtfilt或者改 FIR 线性相位结构。这一步能直接暴露频谱看着对、时域已经错位的隐性 bug。本文还有配套的精品资源点击获取 MIT-BIH心电信号MATLAB插值、滤波与Simulink陷波
简介《课题二 基于MATLAB平台的心电信号分析系统设计及仿真》是一份面向电子信息、生物医学工程等专业课程设计与毕业设计的PDF技术文档适合正在学习数字信号处理、MATLAB编程及Simulink仿真的学生与工程人员参考。文档围绕MIT-BIH心电数据库展开完整设计流程从txt格式心电数据的读取与实际波形还原到线性插值处理使非均匀采样数据等距化再依据心电信号频域特点设计低通、高通及带通滤波器并结合巴特沃斯或切比雪夫滤波器分析幅频特性、零极点分布与阶跃响应。文中还给出50Hz工频陷波器的设计思路以及直接编程静态仿真与Simulink动态建模两条实现路径最后对处理前后的信号做频谱对比分析并整理课程设计报告撰写要求与参考文献。资源包内仅含1个PDF文件大小约1019KB便于在电脑或移动端直接查阅。目前已有357人学习适合希望系统掌握心电信号预处理、滤波器设计与仿真验证方法的读者按步骤实践。1. 心电信号不是丢进滤波器就能出结果MIT-BIH 的 txt 文件长得挺朴素第一列采样时刻第二列 MLII 导联幅值第三列 V5 导联幅值Tab 分隔。坑不在读取而在时间轴——它不是等间距的采集丢点、标注跳变都会让 Δt 在 0.8ms 到 3ms 之间漂移。而 butter、cheby1、ellip 这些函数默认输入是按 Fs 均匀采样的序列你拿非均匀数据直接滤通带边缘会被抹平ST 段那几十微伏的细节直接淹没在数值误差里。所以这个课题真正的执行顺序是读数据并跳过说明行 → 一次线性插值成 1ms 等间隔 → 设计模拟原型滤波器 → 双线性变换离散化 → 50Hz 陷波 → Simulink 里跑动态仿真 → 用 FFT 前后对比频谱反推参数。适合做信号与系统课程设计、DSP 大作业以及刚上手 MATLAB 信号处理工具箱的人。2. MIT-BIH txt 读取与 1ms 等间隔线性插值2.1 三列数据读取与说明行跳过原始文件前两行是文字说明不是数值。要求是严禁手动删除必须用程序忽略。老教材常用textread但它在新版本里已经不推荐我一般直接用textscan配fgetl。fid fopen(100.txt, r); for k 1:2 fgetl(fid); % 只读走前两行说明文字不解析 end C textscan(fid, %f %f %f); % 三列时间、MLII、V5 fclose(fid); t C{1}; % 采样时刻单位秒 x C{2}; % MLII 导联幅值 % C{3} 是 V5同一心拍的另一导联分析只留一路即可textscan返回的是 cell 数组所以要用C{1}取列。如果你把%f %f %f写成%f%f%f遇到 Tab 分隔同样能解析但列数对不上时会静默错位建议保留空格。读完后先做一次体检any(isnan(x))如果返回 1基本就是说明行没跳干净或者文件里混了空行。常见做法是再补一句x x(~isnan(x)); t t(~isnan(t));兜底。2.2 一次线性插值的两种实现插值公式不复杂相邻两点时间差 Δt需要的插值点数 N Δt/0.001幅值步进 ΔA/N然后从 tj-1 ti-1、Aj-1 Ai-1 逐点推。手写循环更能体现原理也方便答辩时讲清楚。Ts 0.001; % 目标采样间隔 1ms tq t(1); xq x(1); % 预置首点 for i 2:length(t) dt t(i) - t(i-1); N round(dt / Ts); % 两点之间需要插几个点 if N 1 % 已经够密不插直接接上 tq(end1,1) t(i); xq(end1,1) x(i); continue; end for j 1:N-1 tq(end1,1) t(i-1) j*Ts; xq(end1,1) x(i-1) (x(i)-x(i-1))*j/N; end tq(end1,1) t(i); xq(end1,1) x(i); endround而不是floor是关键用floor会让累计时间轴越走越慢10 秒数据末尾能偏出几十毫秒。上面这种end1追加写法在大数据量下效率很差正式跑之前把tq、xq预分配成zeros(sum(round(diff(t)/Ts))1, 1)会快一个量级。工程里更省事的是直接调interp1tq t(1):Ts:t(end); xq interp1(t, x, tq, linear); % 默认就是 linear写出来更清楚两者结果在数值上一致差别只是interp1对t中重复时刻点会报警告而手写循环会算出 N0 然后跳过。如果数据里确实有时间戳重复先[t, idx] unique(t); x x(idx);去重再插。2.3 插值结果的校验插完不要急着画图交差。先验三件事max(diff(tq))应该等于 0.001浮点误差范围内如果出现 0.002说明有段原始数据 Δt 小于 0.5ms 被round吞掉了。length(tq)约等于t(end)/0.0013600 点原始数据大约插到 10000 点上下。corr(x, interp1(tq, xq, t))接近 1说明插值没有引入明显畸变。fprintf(插值后点数 %d最大间隔 %.6f s\n, numel(tq), max(diff(tq))); plot(t, x, k., tq, xq, r-); % 原始点叠加插值曲线 xlabel(时间 / s); ylabel(幅值 / mV); legend(原始,插值后);把插值后的[tq, xq]另存为ecg_interp.txt用dlmwrite(ecg_interp.txt, [tq xq], delimiter, \t, precision, %.6f)后面 Simulink 和滤波程序都从这里读避免每次重复插值。3. 巴特沃斯与切比雪夫原型及双线性变换3.1 从心电频谱倒推通带阻带指标正常人心电能量集中在 0.05Hz 到 100HzQRS 波主频在 1025HzP 波和 T 波更低。幅度上成人约 1mV胎儿只有 10μV 量级。这意味着两件事低频端要压住呼吸引起的 0.150.5Hz 基线漂移高频端要挡掉肌电和工频。常见指标组合是低通 Fpass100Hz、Fstop102Hz高通 Fpass1Hz、Fstop0.5Hz阻带衰减 40dB 以上。MIT-BIH 的采样率是 360Hz奈奎斯特 180Hz100Hz 通带留了足够过渡带。滤波器类型通带/阻带阶数作用高通巴特沃斯 IIRFpass1Hz, Fstop0.5Hz4去基线漂移低通等波纹 FIRFpass100Hz, Fstop102Hz最小阶去高频肌电陷波巴特沃斯带阻4853Hz4去 50Hz 工频选巴特沃斯是因为通带最平坦心电的 ST 段对通带波动很敏感1dB 的起伏就能让抬高/压低判断出错。切比雪夫 I 型过渡带更陡同样阶数下能压得更狠代价是通带等波纹用在陷波这种窄带场景更划算。3.2 模拟原型与 freqs 幅频特性先在模拟域把原型搭出来画freqs曲线验证指标这一步是答辩里最容易加分的地方。Rp 1; Rs 40; % 通带 1dB阻带 40dB Wp 2*pi*100; Ws 2*pi*130; % 模拟角频率 rad/s [N, Wc] buttord(Wp, Ws, Rp, Rs, s); [bs, as] butter(N, Wc, s); % 模拟低通原型系数 w logspace(0, 3, 800); % 1 rad/s 到 1000 rad/s H freqs(bs, as, w); semilogx(w/(2*pi), 20*log10(abs(H))); grid on; xlabel(频率 / Hz); ylabel(幅度 / dB);buttord里带s才返回模拟域阶数和截止频率不带就是数字域。Wc通常落在Wp和Ws之间别拿它当 Fpass 用。切比雪夫版本只需换成[N, Wc] cheb1ord(Wp, Ws, Rp, Rs, s)和cheby1(N, Rp, Wc, s)注意cheby1要显式传通带波纹 Rp。freqs给的是模拟频率响应看的是原型形状对不对。零极点分布用[z, p, k] buttap(N)拿归一化原型再把p乘Wc去归一化。所有极点都应该落在 s 左半平面只要有一个跑到右半平面这个滤波器就是发散的不稳定系统必须回头查阶数或参数。3.3 双线性变换的预畸变模拟原型到数字滤波器有两条路。脉冲响应不变法时域逼近好但频谱会混叠只适合低通带通双线性变换把整个 s 平面压进单位圆没有混叠代价是频率轴非线性压缩Ω (2/T)·tan(ω/2)。所以数字截止频率必须先做预畸变再传给模拟原型。Fs 360; % MIT-BIH 采样率 fp 100; fsb 130; % 数字域 Hz Wp 2*Fs*tan(pi*fp/Fs); % 预畸变到模拟域 Ws 2*Fs*tan(pi*fsb/Fs); [N, Wc] buttord(Wp, Ws, Rp, Rs, s); [bs, as] butter(N, Wc, s); [bz, az] bilinear(bs, as, Fs); % 双线性离散化不预畸变直接bilinear的后果很隐蔽截止频率会往低频偏100Hz 的设计点实际可能落在 92Hz 附近通带边缘的心电高频成分被削掉QRS 幅度看起来变矮了。检查方法是freqz(bz, az, 1024, Fs)画数字幅频看 -3dB 点是不是还在 100Hz 附近。等效的写法是直接在数字域用buttord(Wp/(Fs/2), Ws/(Fs/2), Rp, Rs)再butterMATLAB 内部帮你做了预畸变结果和上面一致。两条路都可以走但报告里写清楚用了哪一种别混着放参数。4. Simulink 级联二阶节建模与 50Hz 陷波4.1 六阶低通拆成三个二阶节六阶 Butterworth 低通直接展开成一个六阶传递函数系数动态范围大定点实现或长时间仿真容易数值不稳定。标准做法是拆成三个二阶节级联A 0.09036; B [1.2686, 1.0106, 0.9044]; % 各节分母一次项系数 C [0.7051, 0.3583, 0.2155]; % 各节分母常数项 % 每节形式 Hk(z) A*(1 2z^-1 z^-2) / (1 - Bk*z^-1 Ck*z^-2)三个节串起来Numerator填[A 2*A A]Denominator填[1 -B(k) C(k)]。为什么分子统一是A(12z⁻¹z⁻²)因为 A 是总增益按节数均分后的结果(1z⁻¹)²保证每个节在 z-1 处有零点级联后在奈奎斯特频率处形成六重零点这正是 Butterworth 低通的典型特征。滤波器级数的顺序也有讲究。把 Q 值最高、极点最靠近单位圆的那一节放在最后能让前级先衰减掉一部分能量降低末级的动态范围压力。用zp2sos自动排序比自己拍脑袋排更靠谱[z, p, k] butter(N, Wn); sos zp2sos(z, p, k, down, up); % 按极点靠近单位圆程度排序4.2 From Workspace 与 FDATool 导入Simulink 里喂数据有两条路。一条是From Workspace模块变量名填simin结构体格式simin.time tq; simin.signals.values xq; simin.signals.dimensions 1;老版本要求是这种结构体新版本也接受timeseries(tq, xq)。如果报维度不匹配八成是values列向量写成了行向量加个(:)拉直。数据源记得设成0:Ts:t(end)的有限时长默认的inf会让仿真跑不完。另一条是Digital Filter Design模块先把 FDATool 里设计好的1.fda、2.fda、3.fda保存好在模块 Main 页的 Filter 栏选 Imported加载对应文件。FDATool 的参数对应关系是低通选 Lowpass FIR Equiripple Minimum Order填 Fpass100、Fstop102高通选 Highpass IIR Butterworth Minimum Order填 Fpass1、Fstop0.5带阻选 Bandstop IIR Butterworth填 Fstop149、Fpass148、Fpass253、Fstop251。注意FDATool 里 Minimum Order 算出来的阶数和你在命令行用 buttord 算出来的可能差 1因为 FDATool 用的是它自己的过渡带定义。报告里以 FDATool 导出的系数为准命令行版本只用来交叉验证。4.3 50Hz 陷波器的窄带挑战工频干扰是以 50Hz 为中心、带宽小于 1Hz 的窄带噪声。窄带带阻滤波器有个天然矛盾过渡带越窄阶数越高极点越贴近单位圆对系数量化越敏感。用 FDATool 设计 4853Hz 带阻但实际仿真发现 50Hz 只压下去 20dB 左右说明阶数不够。% 直接算 IIR 陷波比 fdatool 更可控 wo 50/(Fs/2); % 归一化中心频率 bw wo/35; % 3dB 带宽约 1.4Hz [b_n, a_n] iirnotch(wo, bw); y filtfilt(b_n, a_n, xq); % 零相位滤波避免 T 波被拖歪iirnotch的第二个参数是归一化带宽wo/35对应大约 1.4Hz改分母就能调窄。带宽小于 1Hz 时幅度响应会出现很深的凹口但相位响应在凹口附近剧烈变化用filter会造成 T 波形态失真所以这里必须用filtfilt双向滤波。代价是输出长度为原始两倍左右的延迟被抵消边界处会有3*(N-1)点的瞬态前后各截掉一段再算频谱。陷波前先确认 50Hz 干扰是不是真的存在。把原始信号做一次 FFT看 50Hz 处有没有明显尖峰。如果尖峰比周围噪声高不到 10dB先别设计陷波器多半是数据本身已经做过工频抑制硬加陷波反而会把 50Hz 附近的心电成分一起削掉。5. 频谱前后对比与参数回调做完滤波最终判断标准只有一个频谱对比。我一般把原始和滤波后的信号各做一次单边幅度谱用对数坐标叠在一起看。function [f, P] spec1(x, Fs) x x(:) - mean(x); % 去直流 N 2^nextpow2(length(x)); % 补到 2 的幂加快 FFT w hann(length(x)); % 加窗抑制泄漏 X fft(x.*w, N); P abs(X)/sum(w); % 幅度归一化 P P(1:N/21); P(2:end-1) 2*P(2:end-1); f (0:N/2)*(Fs/N); end [f1, P1] spec1(xq, 360); [f2, P2] spec1(y, 360); semilogy(f1, P1, k, f2, P2, r); xlabel(频率 / Hz); ylabel(幅度); legend(滤波前,滤波后);看三个位置就知道参数对不对。0.3Hz 附近滤波后应该下降 20dB 以上基线漂移被压住。50Hz凹口深度如果不到 30dB把iirnotch的带宽再收窄一档。100Hz 以上应该快速滚降如果还在缓降说明低通阶数不够。如果频谱看着没什么变化先别急着换滤波器类型按这个顺序回查时间轴单位。MIT-BIH 的 txt 里第一列单位是秒还是毫秒读错了 Fs 就是 0.36Hz频谱整个挤在左端。filtfilt的边界瞬态。开头和结尾各 200 点常常是数值爆炸区算频谱前先y y(200:end-200);。滤波器系数是否用了归一化频率。butter(N, 100)里的 100 是归一化到奈奎斯特的值必须是小于 1 的数写 100 会直接报错或返回全通。零极点位置。zplane(b, a)看极点有没有跑到单位圆外或者贴到 0.999 这种边缘位置后者说明阶数虚高实际已经临界稳定。最后给一个验证闭环的小技巧把滤波前后的信号各自做一次 QRS 检测哪怕只是最简单的阈值法比较 R 峰位置偏移。如果滤波后 R 峰整体平移超过 5ms说明滤波器群延迟没有被补偿重新用filtfilt或者改 FIR 线性相位结构。这一步能直接暴露频谱看着对、时域已经错位的隐性 bug。本文还有配套的精品资源点击获取