心电图脑电图中场电位噪声干扰如何消除 生物医学信号滤波技术提升心电监测精度与脑电分析准确性

为什么噪声总是来捣乱

想象一下，你正在听一场交响乐演奏，突然有人在你耳边大声说话，还开着电视播放广告——你还能听清音乐吗？心电图（ECG）和脑电图（EEG）也是同样的处境。我们体内的心脏电信号和大脑神经元的电化学活动，本身就很微弱，常常只有几十微伏到几百微伏的强度。偏偏外界环境里到处都在发射电磁波，手机、微波炉、日光灯、医院里的大功率设备，甚至医生手里的电刀，都会在我们的测量信号里留下”痕迹”。

更麻烦的是，噪声的频谱经常和信号本身重叠在一起，想靠一个简单的”切掉这个频率”来解决问题，根本行不通。这时候，我们就需要一些更精细、更聪明的手段，把那些不请自来的噪声请出去，同时保护好真正的信号。

场电位噪声从哪来

先搞清楚敌人是谁，才能制定战术。心电和脑电信号中的噪声大致可以分为几类：

工频干扰：这是最常见的一位。国内市电频率是50Hz，欧美是60Hz，这种干扰会在信号里形成一条尖锐的基线，频率固定、幅度稳定，就像是背景里一直嗡嗡作响的空调外机。

肌电噪声：病人紧张、肌肉紧绷、甚至只是轻微颤抖，都会在信号里引入高频成分，通常集中在20Hz到500Hz之间。这类噪声和心电的高频成分重叠，非常麻烦。

基线漂移：病人呼吸、身体移动、电极接触不良，都会导致信号的基线上下浮动。这种漂移频率很低，可能只有0.1Hz到几Hz，但在看P波、T波或者脑电的慢波活动时，会严重影响判读。

电极伪影：电极和皮肤接触不稳定时，会产生一些突然的尖峰或者不规则的波动，看起来很像真实的心电异常，容易误导诊断。

设备噪声：放大器本身、模数转换过程、电源纹波等，也会引入一定水平的本底噪声。

下面这张表格帮你快速梳理一下各类噪声的特点，方便后续选择合适的滤波策略：

噪声类型	频率范围	典型来源	处理难度
工频干扰	50/60Hz及谐波	电源电磁场	低
肌电噪声	20-500Hz	肌肉收缩	中高
基线漂移	0.1-5Hz	呼吸、运动	中
电极伪影	不规则	接触不良	高
设备噪声	全频段	电子元件	中

从经典到现代：滤波技术的演进路线

滤波这个话题，说大也大，说小也小。它本质上就是对信号做”去粗取精”的工作。我们一步一步来，从最基础的开始讲。

陷波滤波器：专治工频干扰

工频干扰是最容易对付的，因为它频率固定、能量集中。陷波滤波器（Notch Filter）就是为它量身定做的。它的原理很简单——在50Hz（或60Hz）的位置挖一个深深的坑，让这个频率的信号几乎完全被衰减掉。

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

# 生成一个简单的测试信号
fs = 500  # 采样率500Hz
t = np.linspace(0, 2, 1000)  # 2秒的信号

# 干净的心电信号模拟（用多个正弦波叠加近似）
clean_ecg = (
    0.5 * np.sin(2 * np.pi * 1.2 * t) +   # P波模拟
    0.8 * np.sin(2 * np.pi * 3.5 * t) +   # QRS波群模拟
    0.3 * np.sin(2 * np.pi * 2.0 * t)     # T波模拟
)

# 叠加50Hz工频干扰
noisy_ecg = clean_ecg + 0.5 * np.sin(2 * np.pi * 50 * t)

# 设计50Hz陷波滤波器
b, a = signal.iirnotch(50.0, 30.0, fs)  # 中心频率50Hz，Q值30
notch_filtered = signal.filtfilt(b, a, noisy_ecg)

# 可视化结果
fig, axes = plt.subplots(3, 1, figsize=(10, 8))

axes[0].plot(t, clean_ecg, 'b-', linewidth=1.5, label='Clean ECG')
axes[0].set_title('Original Clean Signal')
axes[0].legend()
axes[0].set_ylabel('Amplitude')

axes[1].plot(t, noisy_ecg, 'r-', linewidth=1.5, label='Noisy ECG (with 50Hz interference)')
axes[1].set_title('Signal with Power Line Interference')
axes[1].legend()
axes[1].set_ylabel('Amplitude')

axes[2].plot(t, notch_filtered, 'g-', linewidth=1.5, label='After Notch Filter')
axes[2].set_title('After 50Hz Notch Filter')
axes[2].legend()
axes[2].set_ylabel('Amplitude')
axes[2].set_xlabel('Time (s)')

plt.tight_layout()
plt.savefig('notch_filter_demo.png', dpi=150)
plt.show()

运行这段代码后，你会看到原本被50Hz正弦波干扰的信号，在陷波滤波之后变得干净了许多。不过要注意，陷波滤波器的Q值（品质因数）不能设得太极端。Q值太高虽然滤波更精准，但会引入明显的相位失真；Q值太低又可能影响到附近的信号成分。一般情况下，Q值在20到50之间是一个比较稳妥的选择。

带通滤波器：守住信号的”活动范围”

心电信号的主要能量集中在0.5Hz到150Hz之间，脑电信号则主要在0.5Hz到100Hz左右（虽然伽马波能到200Hz，但常规分析一般用到100Hz就够）。带通滤波器的作用就是保留这个范围内的成分，把低于和高于这个范围的噪声一起滤掉。

# 设计带通滤波器
def design_bandpass_filter(lowcut, highcut, fs, order=4):
    """
    设计带通滤波器
    
    参数:
    lowcut: 低频截止频率
    highcut: 高频截止频率
    fs: 采样率
    order: 滤波器阶数
    """
    nyquist = fs / 2
    low = lowcut / nyquist
    high = highcut / nyquist
    b, a = signal.butter(order, [low, high], btype='band')
    return b, a

# 测试不同的带通设置
fs = 500
lowcut = 0.5
highcut = 45

b, a = design_bandpass_filter(lowcut, highcut, fs, order=4)
filtered_ecg = signal.filtfilt(b, a, noisy_ecg)

# 观察频率响应
w, h = signal.freqz(b, a, worN=2000)
freqs = w * fs / (2 * np.pi)

fig, axes = plt.subplots(2, 1, figsize=(10, 8))
axes[0].plot(freqs, 20 * np.log10(np.abs(h)), 'b-')
axes[0].axvline(lowcut, color='r', linestyle='--', label=f'Low cutoff: {lowcut}Hz')
axes[0].axvline(highcut, color='r', linestyle='--', label=f'High cutoff: {highcut}Hz')
axes[0].set_xlabel('Frequency (Hz)')
axes[0].set_ylabel('Magnitude (dB)')
axes[0].set_title('Bandpass Filter Frequency Response')
axes[0].legend()
axes[0].grid(True)

axes[1].plot(t, noisy_ecg, 'r-', linewidth=1, alpha=0.7, label='Noisy')
axes[1].plot(t, filtered_ecg, 'g-', linewidth=1.5, label='Filtered')
axes[1].set_xlabel('Time (s)')
axes[1].set_ylabel('Amplitude')
axes[1].set_title('Bandpass Filtering Result')
axes[1].legend()
axes[1].grid(True)

plt.tight_layout()
plt.savefig('bandpass_filter_demo.png', dpi=150)
plt.show()

这里有一个小细节值得注意：我们在代码里用了 filtfilt 而不是 filter。filtfilt 做的是双向滤波，它先把信号正向滤波一遍，再把反向滤波一遍。这样做的好处是不引入相位延迟——在医疗信号处理中，这一点至关重要。想象一下，如果滤波让QRS波群延迟了几十毫秒，医生在判断心室除极时间时就会出错。

自适应滤波：对付”狡猾”的噪声

有些噪声是固定的，比如工频干扰，用陷波滤波器就能搞定。但有些噪声是”活”的，它会随着病人的活动而变化。肌电噪声就是一个典型例子——病人紧张时幅度变大，放松时变小，频率成分也会变化。这种时候，自适应滤波器就派上用场了。

自适应滤波的核心思想是：有一个参考通道，里面只含有噪声，不含信号。我们用这个参考噪声来动态调整滤波器的参数，尽量把主通道里的噪声减掉。

from scipy import linalg

def adaptive_filter_lms(reference_noise, primary_signal, fs, mu=0.001, filter_length=32):
    """
    使用LMS算法的自适应滤波器
    
    参数:
    reference_noise: 参考噪声信号（只含噪声）
    primary_signal: 主信号（含信号+噪声）
    fs: 采样率
    mu: 步长参数（控制收敛速度）
    filter_length: 滤波器长度
    
    返回:
    filtered_signal: 滤波后的信号
    error: 误差信号
    """
    n_samples = len(primary_signal)
    filtered = np.zeros(n_samples)
    error = np.zeros(n_samples)
    weights = np.zeros(filter_length)
    
    for i in range(filter_length, n_samples):
        # 当前时刻的参考噪声窗口
        x = reference_noise[i-filter_length:i+1][::-1]
        
        # 滤波器的输出估计
        y = np.dot(weights, x)
        
        # 计算误差
        error[i] = primary_signal[i] - y
        
        # 更新权重（LMS算法）
        weights += 2 * mu * error[i] * x
        
        # 滤波后的输出
        filtered[i] = error[i]
    
    return filtered, error

# 模拟场景：主信号是心电+肌电干扰，参考信号是肌电干扰
np.random.seed(42)
t = np.linspace(0, 5, 2500)
fs = 500

# 干净心电信号
clean_ecg = np.zeros_like(t)
for i in range(len(t)):
    phase = (i % 500) / 500  # 心率100bpm
    if phase < 0.1:
        clean_ecg[i] = 0.1 * np.sin(2 * np.pi * phase * 10)
    elif phase < 0.15:
        clean_ecg[i] = 1.0 * np.sin(2 * np.pi * (phase - 0.1) * 20)
    elif phase < 0.2:
        clean_ecg[i] = -0.3 * np.sin(2 * np.pi * (phase - 0.15) * 15)
    elif phase < 0.3:
        clean_ecg[i] = 0.15 * np.sin(2 * np.pi * (phase - 0.2) * 8)
    elif phase < 0.45:
        clean_ecg[i] = 0.25 * np.sin(2 * np.pi * (phase - 0.3) * 5)

# 肌电噪声（参考通道）
emg_noise = np.random.randn(len(t)) * 0.3 * (1 + 0.5 * np.sin(2 * np.pi * 2 * t))

# 主通道信号
primary_signal = clean_ecg + emg_noise

# 自适应滤波
adaptive_result, err = adaptive_filter_lms(emg_noise, primary_signal, fs)

fig, axes = plt.subplots(4, 1, figsize=(12, 10))

axes[0].plot(t, clean_ecg, 'b-', linewidth=1.2)
axes[0].set_title('Clean ECG Signal')
axes[0].set_ylabel('Amplitude')
axes[0].grid(True)

axes[1].plot(t, emg_noise, 'r-', linewidth=1, alpha=0.7)
axes[1].set_title('EMG Noise (Reference Channel)')
axes[1].set_ylabel('Amplitude')
axes[1].grid(True)

axes[2].plot(t, primary_signal, 'g-', linewidth=1)
axes[2].set_title('Primary Signal (ECG + EMG)')
axes[2].set_ylabel('Amplitude')
axes[2].grid(True)

axes[3].plot(t, adaptive_result, 'm-', linewidth=1.2)
axes[3].set_title('Adaptive Filtered Signal')
axes[3].set_ylabel('Amplitude')
axes[3].set_xlabel('Time (s)')
axes[3].grid(True)

plt.tight_layout()
plt.savefig('adaptive_filter_demo.png', dpi=150)
plt.show()

自适应滤波器的效果如何，很大程度上取决于参考通道的质量。如果参考通道里的噪声和主通道里的噪声是完全一致的，那么效果会非常好；如果两者之间有一些差异，那就需要一些额外的处理，比如对齐、缩放等。在临床实际中，有时候会用两个距离较远的电极来构造参考通道，因为那两个电极主要收到的是噪声，信号成分很少。

小波变换：时频分析的高手

前面讲的滤波器都是在频率域做工作的，它们假设信号的频率成分是相对稳定的。但有些噪声是突然出现的、持续时间很短的，比如电极接触不良产生的尖峰伪影。这种噪声用传统滤波器很难处理，因为它们在频域上和信号重叠。

小波变换提供了一个不同的视角。它可以在时间-频率平面上同时分析信号，找出哪些是真正的信号特征，哪些是短暂的噪声尖峰。

import pywt

# 用一个小波来演示去噪
t = np.linspace(0, 1, 500)
fs = 500

# 干净的脑电信号模拟（Alpha波 8-13Hz）
clean_eeg = 0.5 * np.sin(2 * np.pi * 10 * t) + 0.3 * np.sin(2 * np.pi * 12 * t)

# 叠加噪声和尖峰伪影
noise = np.random.randn(len(t)) * 0.1
spike = np.zeros(len(t))
spike[150:160] = np.linspace(0, 2, 10) * np.exp(-((np.arange(10) - 5)**2) / 2)
spike[350:365] = np.linspace(0, -1.5, 15) * np.exp(-((np.arange(15) - 7)**2) / 3)

noisy_eeg = clean_eeg + noise + spike

# 小波去噪
wavelet = 'db4'
level = 5
coeffs = pywt.wavedec(noisy_eeg, wavelet, level=level)

# 阈值去噪
threshold = np.std(coeffs[-1]) * np.sqrt(2 * np.log(len(noisy_eeg)))
denoised_coeffs = [pywt.threshold(c, threshold * 0.5, mode='soft') for c in coeffs]
denoised_eeg = pywt.waverec(denoised_coeffs, wavelet)

fig, axes = plt.subplots(3, 1, figsize=(12, 8))

axes[0].plot(t, clean_eeg, 'b-', linewidth=1.2, label='Clean EEG')
axes[0].set_title('Clean EEG Signal (Alpha Waves)')
axes[0].legend()
axes[0].grid(True)

axes[1].plot(t, noisy_eeg, 'r-', linewidth=1, label='Noisy EEG')
axes[1].set_title('Noisy EEG with Spikes')
axes[1].legend()
axes[1].grid(True)

axes[2].plot(t, denoised_eeg, 'g-', linewidth=1.2, label='Denoised EEG')
axes[2].set_title('Wavelet Denoising Result')
axes[2].legend()
axes[2].set_xlabel('Time (s)')
axes[2].grid(True)

plt.tight_layout()
plt.savefig('wavelet_denoising.png', dpi=150)
plt.show()

小波变换的好处是它不是”一刀切”地滤掉某个频段，而是可以根据信号在不同尺度上的表现来智能判断。对于脑电信号来说，这种时频局部化的能力特别有价值，因为脑电本身就有很多瞬态成分，比如.event-related potentials（事件相关电位），如果用一个宽频的带通滤波器，可能会把这些重要的特征模糊掉。

盲源分离：把混合信号”拆开”

有时候我们面对的情况更复杂：不仅有心电和脑电，还混进了肌电、基线漂移等各种噪声，而且我们没有参考通道。这种情况下，盲源分离（Blind Source Separation, BSS）技术可以大显身手。

独立成分分析（ICA）是最常用的盲源分离方法之一。它的假设是观测到的信号是若干个独立源信号的线性混合，我们需要找到一个分离矩阵，把这些源信号拆开。

from sklearn.decomposition import FastICA
import numpy as np
import matplotlib.pyplot as plt

# 模拟混合信号场景
fs = 250
t = np.linspace(0, 4, int(fs * 4))

# 源信号
s1 = np.sin(2 * np.pi * 1.2 * t) * np.sin(2 * np.pi * 0.3 * t)  # 低频成分
s2 = np.random.randn(len(t)) * 0.5  # 噪声
s3 = np.sin(2 * np.pi * 8 * t) + 0.5 * np.sin(2 * np.pi * 12 * t)  # 脑电Alpha成分

S = np.vstack([s1, s2, s3])

# 混合矩阵（模拟电极接收到的混合信号）
A = np.array([[0.8, 0.3, 0.2],
              [0.5, 0.6, 0.1],
              [0.2, 0.4, 0.9]])

# 观测信号
X = A @ S

# 使用ICA分离源信号
ica = FastICA(n_components=3, random_state=42)
S_estimated = ica.fit_transform(X.T).T

fig, axes = plt.subplots(4, 1, figsize=(12, 10))

axes[0].plot(t, S[0], 'b-', linewidth=1)
axes[0].set_title('Source 1: Low Frequency Component')
axes[0].set_ylabel('Amplitude')

axes[1].plot(t, S[1], 'r-', linewidth=1)
axes[1].set_title('Source 2: Random Noise')
axes[1].set_ylabel('Amplitude')

axes[2].plot(t, S[2], 'g-', linewidth=1)
axes[2].set_title('Source 3: Brain Alpha Waves')
axes[2].set_ylabel('Amplitude')

axes[3].plot(t, S_estimated[2], 'm-', linewidth=1)
axes[3].set_title('ICA Estimated Source 3')
axes[3].set_ylabel('Amplitude')
axes[3].set_xlabel('Time (s)')

plt.tight_layout()
plt.savefig('ica_demo.png', dpi=150)
plt.show()

在脑电分析中，ICA经常被用来分离眼电（EOG）和肌电（EMG）伪影。眼电伪影的特点是频率很低（0.5-4Hz），而且在特定的导联上幅度很大。通过ICA分解后，我们可以识别出哪些独立成分主要来自眼睛活动，然后把这些成分从信号中去除，保留真正的大脑皮层活动。

组合拳：多种滤波技术的协同使用

在实际的临床应用或科研中，单一滤波技术往往不够用。一个典型的信号处理流程，会把多种技术组合起来，形成一道”过滤防线”。

下面这个流程图展示了从原始信号到最终可用信号的完整处理链路：

原始信号
    │
    ▼
┌─────────────────────┐
│  第一步：预处理       │
│  • 检查信号质量     │
│  • 标记严重伪影     │
│  • 剔除坏导联       │
└──────────┬──────────┘
           │
           ▼
┌─────────────────────┐
│  第二步：去除高频噪声│
│  • 低通滤波 (100-150Hz) │
│  • 目的：去除肌电高频部分  │
└──────────┬──────────┘
           │
           ▼
┌─────────────────────┐
│  第三步：去除工频干扰  │
│  • 陷波滤波 (50/60Hz) │
│  • 或使用频谱抠除     │
└──────────┬──────────┘
           │
           ▼
┌─────────────────────┐
│  第四步：去除基线漂移  │
│  • 高通滤波 (0.5-1Hz) │
│  • 或使用样条拟合     │
└──────────┬──────────┘
           │
           ▼
┌─────────────────────┐
│  第五步：高级处理     │
│  • ICA去除眼电/肌电  │
│  • 小波去噪          │
│  • 自适应滤波        │
└──────────┬──────────┘
           │
           ▼
        可用信号

每一步都有它的必要性。低通滤波在前，是因为高频肌电噪声的能量很大，如果不先压一压，后面的处理可能会被它干扰。陷波滤波放在中间，是因为它的相位特性比较特殊，放太前或太后都可能影响整体效果。高通滤波去除基线漂移放在后面，是因为基线漂移的频率很低，如果在第一步就高通滤波，可能会把一些低频的信号特征也滤掉。

实战：用Python构建完整的心电滤波流程

下面这个例子展示了一个相对完整的心电滤波流程，你可以直接运行看看效果：

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
import pywt

class ECGFilter:
    """
    心电图预处理滤波器类
    
    功能:
    1. 工频干扰消除
    2. 基线漂移去除
    3. 肌电噪声抑制
    4. 小波去噪
    """
    
    def __init__(self, fs, notch_freq=50):
        """
        初始化滤波器
        
        参数:
        fs: 采样率
        notch_freq: 工频频率 (50Hz或60Hz)
        """
        self.fs = fs
        self.notch_freq = notch_freq
        
    def remove_powerline(self, signal_data):
        """
        去除工频干扰
        使用陷波滤波器
        """
        b, a = signal.iirnotch(self.notch_freq, 30, self.fs)
        return signal.filtfilt(b, a, signal_data)
    
    def remove_baseline_drift(self, signal_data):
        """
        去除基线漂移
        使用高通滤波 + 样条拟合
        """
        # 方法一：高通滤波
        lowcut = 0.5
        nyquist = self.fs / 2
        low = lowcut / nyquist
        b, a = signal.butter(2, low, btype='high')
        baseline = signal.filtfilt(b, a, signal_data)
        
        # 方法二：样条拟合基线
        # 找到信号的极小值和极大值点
        from scipy.interpolate import UnivariateSpline
        
        # 简单方法：移动窗口取均值
        window_size = int(self.fs * 0.5)  # 0.5秒窗口
        baseline_spline = self._spline_baseline(signal_data, window_size)
        
        return signal_data - baseline_spline
    
    def _spline_baseline(self, signal_data, window_size):
        """使用样条拟合基线漂移"""
        # 找到局部极小值点
        from scipy.signal import find_peaks
        
        # 反转信号找极大值
        inverted = -signal_data
        peaks, _ = find_peaks(inverted, distance=window_size//2)
        
        if len(peaks) < 2:
            return np.zeros_like(signal_data)
        
        # 样条拟合
        t = np.arange(len(signal_data))
        spline = UnivariateSpline(t[peaks], signal_data[peaks], s=len(signal_data)*0.1)
        return spline(t)
    
    def remove_emg_noise(self, signal_data):
        """
        抑制肌电噪声
        使用低通滤波 + 小波去噪
        """
        # 低通滤波去除高频肌电成分
        highcut = 40  # 心电分析通常40Hz就够了
        nyquist = self.fs / 2
        high = highcut / nyquist
        b, a = signal.butter(4, high, btype='low')
        filtered = signal.filtfilt(b, a, signal_data)
        
        # 小波去噪进一步处理
        wavelet = 'db4'
        level = 4
        coeffs = pywt.wavedec(filtered, wavelet, level=level)
        
        # 设定阈值
        threshold = np.std(coeffs[-1]) * np.sqrt(2 * np.log(len(filtered)))
        
        # 只对小波系数阈值处理，保留近似系数
        denoised_coeffs = [coeffs[0]]  # 保留最低频的近似系数
        for c in coeffs[1:]:
            denoised_coeffs.append(pywt.threshold(c, threshold * 0.5, mode='soft'))
        
        return pywt.waverec(denoised_coeffs, wavelet)
    
    def process(self, signal_data):
        """
        完整的滤波流程
        """
        print("开始心电图信号处理...")
        
        # 第一步：去除工频干扰
        print("  - 去除工频干扰...")
        result = self.remove_powerline(signal_data)
        
        # 第二步：去除基线漂移
        print("  - 去除基线漂移...")
        result = self.remove_baseline_drift(result)
        
        # 第三步：抑制肌电噪声
        print("  - 抑制肌电噪声...")
        result = self.remove_emg_noise(result)
        
        print("处理完成!")
        return result


def generate_sample_ecg(fs=500, duration=5):
    """生成模拟的心电信号"""
    t = np.linspace(0, duration, int(fs * duration))
    
    ecg = np.zeros_like(t)
    heart_rate = 75  # bpm
    beat_duration = 60 / heart_rate
    
    for i in range(len(t)):
        phase = (i / fs) % beat_duration
        
        # P波
        if 0 < phase < 0.08:
            ecg[i] = 0.15 * np.sin(np.pi * phase / 0.08)
        # PR段
        elif 0.08 <= phase < 0.12:
            ecg[i] = 0
        # Q波
        elif 0.12 <= phase < 0.14:
            ecg[i] = -0.1 * np.sin(np.pi * (phase - 0.12) / 0.02)
        # R波
        elif 0.14 <= phase < 0.18:
            ecg[i] = 1.0 * np.sin(np.pi * (phase - 0.14) / 0.04)
        # S波
        elif 0.18 <= phase < 0.22:
            ecg[i] = -0.2 * np.sin(np.pi * (phase - 0.18) / 0.04)
        # ST段
        elif 0.22 <= phase < 0.32:
            ecg[i] = 0.02
        # T波
        elif 0.32 <= phase < 0.45:
            ecg[i] = 0.25 * np.sin(np.pi * (phase - 0.32) / 0.13)
        # TP段
        else:
            ecg[i] = 0
    
    return t, ecg


# 运行演示
fs = 500
t, clean_ecg = generate_sample_ecg(fs)

# 叠加噪声
noisy_ecg = clean_ecg + 0.3 * np.sin(2 * np.pi * 50 * t) + \
            np.random.randn(len(t)) * 0.1 + \
            0.2 * np.sin(2 * np.pi * 0.3 * t)  # 基线漂移

# 应用滤波器
filter_obj = ECGFilter(fs, notch_freq=50)
filtered_ecg = filter_obj.process(noisy_ecg)

# 可视化
fig, axes = plt.subplots(3, 1, figsize=(12, 9))

axes[0].plot(t, clean_ecg, 'b-', linewidth=1.2)
axes[0].set_title('Clean ECG Signal')
axes[0].set_ylabel('Amplitude (mV)')
axes[0].grid(True)

axes[1].plot(t, noisy_ecg, 'r-', linewidth=1)
axes[1].set_title('Noisy ECG Signal')
axes[1].set_ylabel('Amplitude (mV)')
axes[1].grid(True)

axes[2].plot(t, filtered_ecg, 'g-', linewidth=1.2)
axes[2].set_title('Filtered ECG Signal')
axes[2].set_ylabel('Amplitude (mV)')
axes[2].set_xlabel('Time (s)')
axes[2].grid(True)

plt.tight_layout()
plt.savefig('ecg_filtering_demo.png', dpi=150)
plt.show()

脑电信号的特殊处理挑战

脑电信号比心电更难处理，因为它的频率范围更宽、幅度更小、而且噪声来源更复杂。下面这张表对比了心电和脑电在滤波需求上的差异：

特性	心电（ECG）	脑电（EEG）
主要频段	0.5-45Hz	0.5-100Hz+
幅度	0.1-5mV	10-100μV
主要干扰	工频、肌电、基线漂移	工频、眼电、肌电、心电
关键特征	P波、QRS波、T波	Alpha、Beta、Theta、Delta波
采样率需求	250-500Hz	256-1000Hz

脑电的一个特殊问题是眼电（EOG）干扰。眼睛眨动、眼球转动会在额叶导联产生很大的伪影，幅度甚至可以超过脑电本身。解决这个问题的方法之一就是使用ICA，把眼电成分分离出来然后去除。

def remove_eog_artifacts(eeg_data, fs, n_components=4):
    """
    使用ICA去除眼电伪影
    
    参数:
    eeg_data: 多导联脑电数据 (n_channels, n_samples)
    fs: 采样率
    n_components: ICA分解的成分数
    
    返回:
    eeg_clean: 去伪影后的脑电
    """
    from sklearn.decomposition import FastICA
    
    # 预处理：标准化
    eeg标准化 = (eeg_data - np.mean(eeg_data, axis=1, keepdims=True)) / np.std(eeg_data, axis=1, keepdims=True)
    
    # ICA分解
    ica = FastICA(n_components=n_components, random_state=42)
    sources = ica.fit_transform(eeg标准化.T).T
    
    # 识别眼电成分（通常在前几个导联幅度最大）
    eog_components = []
    for i in range(n_components):
        # 检查每个成分在各个导联的方差
        variance_profile = np.var(sources[i], axis=1)
        # 如果额叶导联（前2-3个）方差明显大于其他导联，可能是眼电
        frontal_var = np.mean(variance_profile[:3])
        total_var = np.mean(variance_profile)
        if frontal_var > 2 * total_var:
            eog_components.append(i)
    
    print(f"检测到{len(eog_components)}个眼电成分: {eog_components}")
    
    # 去除眼电成分
    sources_clean = sources.copy()
    for comp_idx in eog_components:
        sources_clean[comp_idx] = 0
    
    # 重建信号
    eeg标准化_clean = ica.inverse_transform(sources_clean.T).T
    
    # 还原到原始尺度
    eeg_clean = eeg标准化_clean * np.std(eeg_data, axis=1, keepdims=True) + np.mean(eeg_data, axis=1, keepdims=True)
    
    return eeg_clean

实际应用中需要注意的细节

理论讲了一大堆，但在实际做项目或者处理真实数据的时候，还有很多细节需要注意。

采样率的选择：采样率不能太低，否则会丢失信号的高频成分，也会引起混叠。心电一般250-500Hz就够了，脑电推荐500Hz以上，如果要分析高伽马波（>80Hz），可能需要1000Hz甚至更高。

滤波器的相位特性：在医疗应用中，相位失真是一个严重的隐患。IIR滤波器（如Butterworth）会有相位延迟，虽然用filtfilt做双向滤波可以消除相位失真，但会增加计算量。FIR滤波器可以设计成线性相位的，但阶数通常很高。根据实际需求权衡取舍。

滤波器的稳定性：尤其是IIR滤波器，要检查极点是否在单位圆内。对于递归滤波器，数值精度也可能导致不稳定，特别是在处理长信号时。

过拟合的风险：滤波太猛会把信号的有用成分也滤掉。比如高通滤波切得太高，可能会把P波的低频部分也滤掉，影响心电的判读。建议滤波前后都对比一下波形，确保没有过度处理。

实时 vs 离线：如果是实时监测，滤波器的延迟必须控制在一个可接受的范围内。IIR滤波器延迟小，但相位失真；FIR滤波器可以做到线性相位，但延迟较大。这是一个需要权衡的问题。

总结与展望

心电和脑电的噪声消除是一个需要综合考量多个因素的问题。从经典的陷波滤波、带通滤波，到自适应滤波、小波变换、ICA，每种技术都有自己的适用场景和优缺点。实际应用中，往往需要组合多种方法，根据具体的噪声特征和信号特点来设计处理流程。

未来的趋势是朝着智能化发展。深度学习模型已经在一些伪影去除任务上展现了不错的效果，比如用卷积神经网络自动识别和去除眼电伪影，或者用生成对抗网络来重建被噪声污染的脑电 segment。不过，传统信号处理方法在可解释性和计算效率上仍有优势，短期内不太可能被完全取代。

最好的策略是：用传统方法做好预处理和基础滤波，然后用数据驱动的方法做精细化处理，两者结合才能达到最佳效果。希望这篇文章能帮你建立起一个比较完整的认知框架，在实际工作中根据具体情况灵活运用。