简介：面向生物医学信号处理与胎儿心电监测方向的学习者，这份MATLAB实现围绕LMS自适应滤波提取FECG展开，可用于理解如何从含母体心电及噪声的腹部观测信号中分离胎儿心电成分。包内共2个文件，以m脚本和dat数据样本为主，压缩包约46KB，体量轻便，适合在MATLAB环境中直接读取数据并复现自适应滤波实验。已有917人学习，说明其在相关课程设计、算法入门与实验验证中具有一定参考价值。通过脚本与实测数据的配合，读者可以观察LMS权值更新、误差收敛和滤波输出变化，进而分析算法参数对胎儿心电增强效果的影响，并在此基础上做预处理、收敛步长调整或结果可视化的二次尝试。整体更偏向基础算法演示与小规模实验样本，适合具备一定数字信号处理基础、希望快速接触FECG提取流程的本科生和初级研究者参考。

1. 从腹壁混合信号里捞出胎儿心电：这套 LMS 自适应滤波 MATLAB 算法解决的具体问题

产科监护里有个老问题：贴在孕妇腹壁上的电极，采到的是母体心电、胎儿心电和肌电噪声的叠加。母体 QRS 幅度通常是胎儿的 5 到 20 倍，两者频谱又大幅重叠，普通带通滤波器切不干净，一刀下去连胎儿成分一起削掉。这套以 LMS 自适应滤波为核心的 MATLAB 算法，思路换了个方向——再取一路胸导作为参考，让自适应滤波器持续估计腹导里的母体心电，再从原始信号里减掉它，残差里剩下的那点微弱波形就是 FECG。它适合三类人：做生物医学信号处理课程设计、要把算法跑出可见波形、想拿一份能跑通的代码当自适应滤波入门模板。

2. 信号链路拆解：参考通道、时延对齐与 LMS 抵消器的三个输入

动手写代码前先把信号模型理清，这一步决定了后面参数怎么调。很多人直接开写循环，结果跑出来一条平线，回头查半天发现是参考信号选错了，或者根本没收敛。先把模型摆正，流程就能顺畅推进。

2.1 主输入、参考输入与残差各自的物理含义

自适应噪声抵消器（ANC）里有三个信号，命名容易混，先定死：

写成公式就是d(n) = s(n) + h(n) * m(n) + v(n)，其中h(n)是从心脏到腹壁这条传播路径的冲激响应，未知、时变、且因人而异。LMS 滤波器的权向量w干的事情就是在线逼近这个h(n)：只要w能把参考通道卷成腹壁上的母体心电形状，减法的结果里母体成分自然被压下去。

这里有个容易忽略的点：h(n)不是简单的比例缩放，它带了时延和相位失真。母体心电从心脏传到胸导几乎没延迟，传到腹壁要经过一段组织，两路信号差个十几到几十毫秒很常见。所以自适应滤波器的阶数L必须足够覆盖这段“信道记忆”，光靠一个系数乘一乘是对消不掉的。我一般按经验取L = 32（采样率 1 kHz 时约合 32 ms 的记忆窗），如果对齐后残差里还有明显的母体 QRS 残留，就往上加。

2.2 时延对齐为什么不能跳过

参考信号和主输入如果没对齐，LMS 要多花好几倍的阶数去“造”出一个延迟，收敛变慢，稳态失调变大，直观表现就是前几秒输出完全没法看。对齐也不难做，用互相关搜索主峰位置就行：

function [ab, th, delay] = align_channels(ab, th, maxLag) % 用互相关估计胸导相对腹导的母体 QRS 时延，并把两路拉齐 ab = ab(:) - mean(ab); % 去直流，避免直流分量主导互相关 th = th(:) - mean(th); [c, lags] = xcorr(ab, th, maxLag); [~, idx] = max(abs(c)); delay = lags(idx); % 峰值位置即两路的相对时延（样本点） if delay > 0 th = [th(delay+1:end); zeros(delay,1)]; % 胸导提前，补偿滞后 elseif delay < 0 th = [zeros(-delay,1); th(1:end+delay)]; end end

maxLag是搜索半径，取值要大于物理上可能的最大传播时延，采样率 1 kHz 时给 200 个点（200 ms）绰绰有余；给太大会让互相关在错误的周期峰上锁定，母体心率 60 bpm 时心动周期是 1000 ms，只要搜索窗远小于一个心动周期就不会串到隔壁峰。函数末尾那段移位逻辑里的符号方向，建议第一次用时先把delay打印出来，拿一路已知的合成信号对一遍，别照抄——这个地方出错，后面所有分析都会白费。

2.3 步长、阶数与收敛性的定量关系

LMS 的权重更新式是w(n+1) = w(n) + 2μ·e(n)·x(n)，μ就是步长。它决定了两个互相矛盾的指标：收敛速度和稳态失调。工程上常用两条经验边界来框定：

• 收敛条件：0 < μ < 1 / (L · σx²)，σx²是参考信号功率。胸导信号归一化到 ±1 之后，L = 32对应μ的上界大约 0.03，实际取 0.001～0.005 之间比较稳。
• 时间常数：τ ≈ 1 / (2μ·λ_avg)，λ_avg是参考信号自相关矩阵的平均特征值。想让它 2 秒内收敛，就把τ对应的样本数换算回去反推μ。

这两条式子不是让你精确计算，而是给你一个调试方向：跑出来发现前 10 秒还在剧烈震荡，说明μ太大或者L太大；跑出来发现权重压根不动，残差和原始腹导几乎一样，那就是μ太小。我通常从μ = 0.002、L = 32起步，先看权能量收敛曲线，再回头调。别一上来就调阶数，阶数一变大收敛条件跟着收紧，两个参数一起动会把自己绕进去。

3. 工程落地：MATLAB 目录结构、主脚本六步流程与参数配置

算法想清楚了，剩下的是把它组织成别人能看懂、自己能反复跑的形式。散装脚本最大的问题是改一次参数就得回头翻半天，所以我习惯把参数、算法、评估拆成三层。

3.1 目录与文件职责划分

fecg_lms/ ├── main_fecg.m % 唯一入口，六步串起来 ├── config.m % 所有可调参数集中在这里 ├── src/ │ ├── preprocess.m % 去直流、去工频、重采样到统一 fs │ ├── align_channels.m % 互相关时延对齐 │ ├── lms_anc.m % LMS 自适应抵消主函数 │ ├── nlms_anc.m % 归一化版本，用于幅度漂移大的记录 │ ├── detect_fqrs.m % 胎儿 QRS 检测（带通+平方+滑窗积分） │ └── eval_snr.m % 母体抑制比、胎儿段信噪比 ├── data/ │ └── record_01.mat % 两路导联记录，含 fs 字段 └── out/ % 图、指标表、中间变量都落这里

这种拆法有个好处：lms_anc.m只依赖两个向量和两个标量，完全不知道数据从哪来，换个数据集不用改一行算法代码。写论文或做对比实验时，直接从src/里取函数当黑盒调用就行。

3.2 主脚本的六步执行流程

main_fecg.m我刻意写得很短，短到一眼能扫完，因为流程写长了以后改一处忘一处，我遇到过太多次。

%% main_fecg.m —— LMS 提取胎儿心电的完整流程 clc; clear; close all; cfg = config(); % 1) 读参数 S = load(fullfile(cfg.dataDir, cfg.recordFile)); [ab, th, fs] = preprocess(S.abdominal, S.thoracic, S.fs, cfg); % 2) 预处理 [ab, th, dly] = align_channels(ab, th, cfg.maxLag); % 3) 时延对齐 fprintf('通道时延 = %d 样本 (%.1f ms)\n', dly, 1000*abs(dly)/fs); [e, w, wnorm] = lms_anc(ab, th, cfg.mu, cfg.L); % 4) LMS 抵消 fqrs = detect_fqrs(e, fs, cfg); % 5) 胎儿 QRS 定位 metrics = eval_snr(e, ab, fs, cfg); % 6) 指标与出图 save(fullfile(cfg.outDir, 'result.mat'), 'e', 'w', 'fqrs', 'metrics');

每一步的输入输出都落到变量上，中间想看哪一步，直接在对应行打断点，不用重跑全流程。第 3 步打印的时延值我建议一直留着，它是判断对齐是否合理的第一个信号，正常应该在 0～50 ms 之间，跑出三四百毫秒基本可以断定互相关锁错周期了。

第 4 步返回三个值：e是残差（也就是提取出的胎儿心电雏形），w是最终权向量，wnorm是逐样本的权能量‖w‖²。第三个值调试时特别有用，它收敛到平台期说明滤波器稳定了，一直往上飘说明步长越界开始发散。

3.3 参数集中管理与结果落盘

config.m是一个返回结构体的函数，不写脚本级变量，避免和主脚本的工作区互相污染：

function cfg = config() cfg.fs = 1000; % 统一采样率 Hz cfg.L = 32; % 自适应滤波器阶数（覆盖信道记忆长度） cfg.mu = 0.002; % LMS 步长，需满足 mu < 1/(L*var(x)) cfg.maxLag = 200; % 互相关时延搜索半径（样本点） cfg.bp = [1 40]; % 后处理带通频带 Hz cfg.win = 0.12; % 滑窗积分宽度 s cfg.notch = 50; % 工频陷波中心 Hz，按现场电源频率改 cfg.dataDir = './data'; cfg.recordFile = 'record_01.mat'; cfg.outDir = './out'; end

最有价值的其实是fs和cfg.notch。国内现场交流电频率按 50 Hz，某些设备导出时已经被硬件滤过一遍，再陷波反而会把胎儿 QRS 的高频成分削掉一块；采样率不一致更常见，不同设备导出的记录从 250 Hz 到 1 kHz 都有，preprocess.m里统一用resample拉到cfg.fs，后面所有以样本点计的参数（阶数、滑窗宽度、最小峰间距）才有意义。这一步不做，换个数据集就得重新调一遍参，很折腾。

4. 核心代码逐段讲：逐样本更新的 LMS 主循环与 NLMS 归一化实现

真正干活的就十几行，但每一行都有讲究。

4.1 LMS 主循环：从参考快照到残差输出

function [e, w, wnorm] = lms_anc(d, x, mu, L) % LMS 自适应噪声抵消：用参考 x 估计并抵消主输入 d 中的母体心电 % d 主输入（腹导），列向量 % x 参考输入（胸导），列向量 % mu 步长 % L 滤波器阶数 % e 残差，近似等于胎儿心电 N = numel(d); d = d(:); x = x(:); w = zeros(L, 1); % 权向量初值置零 e = zeros(N, 1); wnorm = zeros(N, 1); xb = [zeros(L-1, 1); x]; % 前补零，保证 n=1 时也有完整快照 for n = 1:N xn = xb(n+L-1 : -1 : n); % 当前参考快照：[x(n); x(n-1); ...; x(n-L+1)] y = w.' * xn; % 母体心电的线性估计 e(n) = d(n) - y; % 残差 = 主输入 - 母体估计 w = w + 2*mu*e(n)*xn; % LMS 权重更新 wnorm(n) = w.' * w; % 记录权能量用于观察收敛 end end

逐行看。xb那次前补零是整个实现里最容易写错的地方：如果不补零，用x(n-L+1:n)取快照在n < L时会越界或者被 MATLAB 自动截断，前 L-1 个样本的输出就变成无效值，图上表现为开头一小段剧烈跳动。补零之后，xb(n+L-1:-1:n)在n = 1时取到[x(1); 0; ...; 0]，在n = 2时取到[x(2); x(1); 0; ...; 0]，索引和物理含义完全对得上，写完之后顺手在n = 1和n = L两处打断点验一遍，比事后查半天强。

y = w.' * xn是滤波器的输出，物理上就是“当前估计出的母体心电在腹壁上的样子”。注意它并不是参考信号本身，而是参考信号经过当前权向量卷积之后的结果，这也是自适应滤波和固定滤波的本质区别。

e(n) = d(n) - y是残差。这里有个反直觉的地方：残差不是误差，它同时是算法的输出。LMS 的更新准则是最小化E[e²]，而e里剩下的主要成分就是胎儿心电加上没被抵消掉的噪声，所以最小化残差功率恰好等价于“把母体成分尽力压掉”。

w = w + 2*mu*e(n)*xn是标准的 LMS 更新式。系数 2 来自对e²求梯度，别随手改成 1，改了以后mu的物理含义跟着变，你按公式算出来的边界就不准了。最后一行记录‖w‖²，它比残差曲线更早暴露发散问题：残差开始异常之前，权能量通常已经在往上飘了。

4.2 NLMS 归一化版本：把步长从手动调参变成可复现参数

固定步长的 LMS 有个明确弱点：腹壁信号幅度会随呼吸和体位漂移，同一段记录里参考信号功率可能差三四倍，mu按前半段调好了，后半段就可能越界。归一化 LMS 把步长按当前快照能量做了缩放：

function [e, w] = nlms_anc(d, x, mu, L, eps_reg) % 归一化 LMS：步长按当前参考快照能量自适应缩放 if nargin < 5, eps_reg = 1e-6; end % 防止除零 N = numel(d); d = d(:); x = x(:); w = zeros(L, 1); e = zeros(N, 1); xb = [zeros(L-1, 1); x]; for n = 1:N xn = xb(n+L-1 : -1 : n); e(n) = d(n) - w.' * xn; w = w + mu * e(n) * xn / (xn.'*xn + eps_reg); % 除以快照能量 end end

关键差别在最后一行：分母xn.'*xn是当前 L 点快照的能量，信号强的时候分母大、等效步长自动变小，信号弱的时候反过来。eps_reg这个正则项必须留，静默期参考信号接近全零，不加上去就是一个除零，输出直接变成 NaN，之后整个数组被污染，图上一片空白，这种问题查起来最费劲。归一化之后mu的取值范围可以放开到 0.05～0.5，且对信号幅度基本不敏感，换数据集时省掉一整轮重调参。

那为什么还保留 LMS 版本？因为做算法对比实验时，它是基准。论文里报改进效果，对照组用标准 LMS 才说得过去，直接上 NLMS 反而不好解释增益来自哪里。

4.3 后处理链：从残差到可定位的胎儿 QRS

抵消后的残差还不能直接拿去数胎心，里面残余的肌电噪声和 T 波干扰会让峰值检测乱跳。常规做法是走一遍带通、平方、滑窗积分：

function fqrs = detect_fqrs(e, fs, cfg) % 胎儿 QRS 检测：带通 -> 平方 -> 滑窗积分 -> 自适应阈值寻峰 [b, a] = butter(3, cfg.bp/(fs/2), 'bandpass'); ef = filtfilt(b, a, e); % 零相位滤波，QRS 位置不偏移 en = ef.^2; % 平方：压小噪声，放大 QRS w = round(cfg.win * fs); ei = filter(ones(w,1)/w, 1, en); % 滑窗积分，把 QRS 包络成单峰 ei = ei / max(ei + eps); thr = 0.35 * mean(ei(ei > 0.2*max(ei))); % 自适应阈值 [pks, loc] = findpeaks(ei, 'MinPeakHeight', thr, ... 'MinPeakDistance', round(0.25*fs)); % 对应心率上限 240 bpm fqrs = loc(:); end

filtfilt而不是filter，是因为前者做正反向两次滤波，相位为零，QRS 峰的位置不会被滤波器群延迟推着走，后面算 RR 间期才准。带通下限 1 Hz 用来压掉残余的基线漂移，上限 40 Hz 是折中：胎儿 QRS 的主要能量集中在 10～35 Hz，上限再高会把肌电噪声一起放进来。滑窗宽度0.12 s对应成年人心率的 QRS 宽度量级，对胎儿偏宽一点没关系，积分窗的作用是把一个 QRS 波群压成一个单峰，窗口太窄会分裂成两个峰。

MinPeakDistance设成 0.25 s 是有依据的：胎儿心率上限大约 240 bpm，对应 RR 间期 0.25 s，小于这个间距的峰一定是噪声或者同一个 QRS 的双峰。阈值那行用mean而不是max的一半，是因为残差里可能有个别幅度特别大的伪迹，用最大值归一化会把阈值抬得过高，导致后面一大段检不出来。

5. 排查与避坑：跑不出胎儿 QRS 的五个典型翻车点

这套流程我在不同记录上反复跑过，出问题的地方高度集中，下面几条都是实打实遇到的。

5.1 数据与预处理层面的三个坑

现象一：残差里的母体 QRS 比原始腹导还大。原因基本是参考通道的时延补偿方向搞反了，或者两路信号长度不一致导致移位后错位。解决：先把align_channels返回的delay打印出来，用一段合成的、已知时延的双通道信号验证移位方向；同时确认两路信号在同一时刻开始采集，长度不等就先截到等长再对齐。

现象二：结果是一条接近零的直线，什么都看不出来。这不是滤波太强，而是滤波器把信号整体消掉了。多半是参考通道和主输入接成了同一路（比如数据文件里两列存重了，或者变量赋值时手滑写成th = ab）。解决：进preprocess前后各算一次两路的相关系数，正常应该很高但不是 1，接近 1 就说明拿错了。另外检查预处理有没有把某一路整体归零，陷波器参数不合法时会出现这种情况。

现象三：波形里有一层均匀的高频毛刺，胎儿 QRS 淹没在里面。原因是工频干扰没滤干净，或者陷波频率设错。设备在不同现场录的数据，工频可能是 50 Hz 也可能是 60 Hz，config里那个值必须按实际记录改。解决：先对原始信号做一次频谱，看 50/60 Hz 处有没有尖峰，有尖峰再开陷波；同时确认陷波器的品质因数别设太高，Q 值过大在中心频率附近会产生很长的振铃，反而污染 QRS。

5.2 参数与收敛层面的两个坑

现象四：开头几秒输出剧烈震荡，之后才慢慢稳定。这是 LMS 的固有收敛过程，不是 bug。但如果震荡超过十几秒还在继续，就是步长越界了。解决：先看wnorm曲线有没有进入平台期，一直单调上升说明在发散，把mu减半再试；也可以直接换nlms_anc，归一化之后对幅度不敏感，这类问题基本消失。判断收敛快慢的实用办法是把mu从 0.0005 扫到 0.01 画一组曲线，肉眼选那个“收敛够快且稳态不抖”的点。

现象五：胎儿 QRS 检测出来的数量明显偏多，心率算出来常年 200 以上。不是算法的问题，是detect_fqrs里阈值定得太松，把残余母体 T 波或者噪声峰检成了胎心。解决：先看残差波形，如果每个母体 QRS 之后仍有一个明显的小峰，说明母体成分没消干净，回头加阶数或者检查对齐；如果波形本身很干净但检测仍多，把MinPeakDistance提到 0.3 s，或者把阈值系数从 0.35 提到 0.5 再跑一遍对比。这两步我都习惯各存一张图，放在out/里，改参数后一对比就知道是滤波变好了还是检测变松了。

提示：每改一次参数就把result.mat落一次盘，命名带上参数值，别用result_final.mat、result_final2.mat这种。回头要复现某张图时，你会庆幸自己当时多花了这几秒。

6. 二次开发改造：换参考构造、换自适应核、补一套可对比的评估

如果只拿这套代码跑一遍出图，价值有限；真正能用起来的地方在于它留了几个干净的替换点。第一个是参考信号的构造。单路胸导做参考，对消母体心电已经够用，但如果腹壁上母体成分的形态和胸导差得多（多胎、胎位异常或者电极贴偏都可能出现），可以扩成多路参考：胸导加它的若干延迟版本拼成一个向量，喂给同一个lms_anc，等价于人为给滤波器更长的记忆，收敛会明显快一些。改法很简单，把x从列向量换成N×K矩阵，主循环里取xn = xb(n+L-1:-1:n, :)，w改成矩阵，其余不动。

第二个是自适应核的替换。lms_anc和nlms_anc的接口完全一致（d, x, mu, L进，e, w出），所以加一个 RLS 版本只需要新写一个同签名函数，主脚本里换个函数名就能对比。三种核的取舍我一般这么把握：

第三个容易被忽略的是评估环节。现在的eval_snr.m只算母体抑制比和胎儿段信噪比，够用但不够有说服力。想验证效果，最直接的办法是构造一段半合成数据：拿干净的母体心电和一段已知位置的单峰脉冲按比例叠加，脉冲位置就是胎儿 QRS 的真值，跑完流程后直接比对检出率与位置误差。这套半合成验证我吃过亏——早期只看图觉得“波形挺像那么回事”，换一批数据后检出率掉到六成，回头才发现是阈值系数和采样率绑死了，改fs就失效。从那以后我每次改完参数，都强制走一遍半合成数据回归，检出率和位置误差两个数都在预期范围内，才敢往真实记录上跑。

真实记录没有真值，但有个粗略的合理性检查：胎儿心率应该落在 110～160 bpm，且逐拍变化平缓，RR 间期序列的标准差一般不超过 30 ms。跑出来心率 90 或者 220，基本都是检测环节而不是滤波环节的问题，别急着回去调mu。

希望这些能帮你少走几步弯路，把第一张能看的胎儿心电波形跑出来。

本文还有配套的精品资源，点击获取