EEG信号处理实战：用Python实现小波变换降噪（附完整代码）

脑电信号（EEG）是神经科学研究的重要数据来源，但原始EEG信号往往包含大量噪声干扰。本文将带你用Python实现基于db4小波的4级分解与重构，解决高频噪声问题，并对比不同小波基的降噪效果。

1. 理解小波变换的核心优势

传统傅里叶变换只能提供信号的全局频率信息，而小波变换同时具备时频局部化特性。这种特性使其特别适合处理非平稳信号如EEG：

• 时域分辨率：高频成分对应更精细的时间分辨率
• 频域分辨率：低频成分对应更精确的频率分辨率

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

# 生成测试信号
fs = 1000  # 采样率
t = np.linspace(0, 1, fs)
signal = np.sin(2*np.pi*10*t) + 0.5*np.sin(2*np.pi*50*t)

# 短时傅里叶变换
f, t_spec, Sxx = spectrogram(signal, fs)
plt.pcolormesh(t_spec, f, 10*np.log10(Sxx))
plt.ylabel('Frequency [Hz]')
plt.xlabel('Time [sec]')
plt.show()

提示：小波变换的时频窗可自适应调整，高频区域时间窗窄，低频区域频率窗窄

2. 数据准备与预处理

使用Bonn大学公开EEG数据集，包含健康人和癫痫患者的脑电记录。我们主要处理以下两类数据：

1. 健康组（Z系列）：睁眼静息状态
2. 癫痫组（S系列）：癫痫发作期

数据预处理流程：

1. 加载原始数据（采样率173.61Hz）
2. 带通滤波（0.5-40Hz）
3. 去除工频干扰（50Hz陷波）

# 数据加载示例
def load_eeg_data(file_path):
    with open(file_path, 'r') as f:
        data = [float(line.strip()) for line in f]
    return np.array(data)

# 带通滤波实现
from scipy.signal import butter, filtfilt

def butter_bandpass(lowcut, highcut, fs, order=5):
    nyq = 0.5 * fs
    low = lowcut / nyq
    high = highcut / nyq
    b, a = butter(order, [low, high], btype='band')
    return b, a

def bandpass_filter(data, lowcut, highcut, fs, order=5):
    b, a = butter_bandpass(lowcut, highcut, fs, order=order)
    return filtfilt(b, a, data)

3. db4小波分解与重构

Daubechies-4（db4）小波因其与EEG信号的相似性被广泛使用。我们实现4级分解：

1. 原始信号 → A1（近似） + D1（细节）
2. A1 → A2 + D2
3. A2 → A3 + D3
4. A3 → A4 + D4

关键参数对比表：

分解层数	频率范围(Hz)	主要成分
D1	43.4-86.8	高频噪声
D2	21.7-43.4	γ波
D3	10.8-21.7	β波
D4	5.4-10.8	α波
A4	0-5.4	δ/θ波

import pywt

def wavelet_decomposition(signal, wavelet='db4', level=4):
    coeffs = pywt.wavedec(signal, wavelet, level=level)
    return coeffs

def wavelet_reconstruction(coeffs, wavelet='db4'):
    return pywt.waverec(coeffs, wavelet)

4. 阈值去噪实战

小波阈值去噪包含三个关键步骤：

1. 系数选择：保留低频近似系数，处理高频细节系数
2. 阈值计算：采用通用阈值（Universal Threshold）
3. 阈值函数：软阈值或硬阈值

Python实现细节：

def wavelet_denoise(signal, wavelet='db4', level=4, mode='soft'):
    # 分解
    coeffs = wavelet_decomposition(signal, wavelet, level)
    
    # 计算阈值
    sigma = np.median(np.abs(coeffs[-1])) / 0.6745
    uthresh = sigma * np.sqrt(2*np.log(len(signal)))
    
    # 应用阈值
    new_coeffs = []
    new_coeffs.append(coeffs[0])  # 保留近似系数
    
    for i in range(1, len(coeffs)):
        if mode == 'soft':
            new_coeffs.append(np.sign(coeffs[i]) * 
                             np.maximum(0, np.abs(coeffs[i]) - uthresh))
        else:  # hard
            new_coeffs.append(coeffs[i] * (np.abs(coeffs[i]) > uthresh))
    
    # 重构
    return wavelet_reconstruction(new_coeffs, wavelet)

注意：软阈值处理会产生更平滑的结果，但可能丢失部分有效信号；硬阈值保留更多细节但可能引入伪影

5. 效果评估与可视化

我们使用信噪比（SNR）和均方根误差（RMSE）量化降噪效果：

def evaluate_denoising(original, noisy, denoised):
    # 计算信噪比
    def snr(signal, noise):
        return 10*np.log10(np.var(signal)/np.var(noise))
    
    # 计算RMSE
    rmse = np.sqrt(np.mean((original - denoised)**2))
    
    print(f"输入SNR: {snr(original, noisy-original):.2f} dB")
    print(f"输出SNR: {snr(original, denoised-original):.2f} dB") 
    print(f"RMSE: {rmse:.4f}")
    
    # 绘制对比图
    plt.figure(figsize=(12, 6))
    plt.subplot(3,1,1)
    plt.plot(original)
    plt.title("Original Signal")
    
    plt.subplot(3,1,2) 
    plt.plot(noisy)
    plt.title("Noisy Signal")
    
    plt.subplot(3,1,3)
    plt.plot(denoised)
    plt.title("Denoised Signal")
    
    plt.tight_layout()
    plt.show()

不同小波基性能对比：

小波类型	计算速度	平滑性	特征保持
db4	中等	好	优
sym4	快	很好	良
coif4	慢	优	优
haar	最快	差	差

6. 癫痫检测应用实例

将小波变换应用于癫痫发作检测，关键步骤包括：

1. 对EEG信号进行5级小波分解
2. 提取各频带能量作为特征
3. 使用SVM分类器检测异常

def extract_wavelet_features(signal, wavelet='db4', level=5):
    coeffs = pywt.wavedec(signal, wavelet, level=level)
    features = []
    
    for i, c in enumerate(coeffs):
        # 计算每个频带的相对能量
        energy = np.sum(c**2)
        features.append(energy)
    
    # 归一化
    total_energy = np.sum(features)
    return [f/total_energy for f in features]

# 示例：从癫痫数据中提取特征
seizure_data = load_eeg_data('S001.txt')
features = extract_wavelet_features(seizure_data)

在实际项目中，这种方法的分类准确率可达85%以上，比传统频域分析方法提高约15%。

7. 完整代码实现

以下是整合后的完整处理流程：

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

class EEGDenoiser:
    def __init__(self, fs=173.61, wavelet='db4', level=4):
        self.fs = fs
        self.wavelet = wavelet
        self.level = level
        
    def bandpass_filter(self, data, lowcut=0.5, highcut=40):
        nyq = 0.5 * self.fs
        low = lowcut / nyq
        high = highcut / nyq
        b, a = butter(5, [low, high], btype='band')
        return filtfilt(b, a, data)
    
    def notch_filter(self, data, freq=50, Q=30):
        nyq = 0.5 * self.fs
        freq = freq / nyq
        b, a = butter(2, [freq-0.01, freq+0.01], btype='bandstop')
        return filtfilt(b, a, data)
    
    def denoise(self, signal, mode='soft'):
        # 预处理
        signal = self.bandpass_filter(signal)
        signal = self.notch_filter(signal)
        
        # 小波分解
        coeffs = pywt.wavedec(signal, self.wavelet, level=self.level)
        
        # 阈值计算
        sigma = np.median(np.abs(coeffs[-1])) / 0.6745
        uthresh = sigma * np.sqrt(2*np.log(len(signal)))
        
        # 阈值处理
        new_coeffs = [coeffs[0]]
        for i in range(1, len(coeffs)):
            if mode == 'soft':
                new_coeffs.append(np.sign(coeffs[i]) * 
                                 np.maximum(0, np.abs(coeffs[i]) - uthresh))
            else:
                new_coeffs.append(coeffs[i] * (np.abs(coeffs[i]) > uthresh))
        
        # 重构
        return pywt.waverec(new_coeffs, self.wavelet)
    
    def plot_decomposition(self, signal):
        coeffs = pywt.wavedec(signal, self.wavelet, level=self.level)
        plt.figure(figsize=(12, 8))
        
        for i, c in enumerate(coeffs):
            plt.subplot(len(coeffs), 1, i+1)
            plt.plot(c)
            plt.title(f'Level {i}' if i == 0 else f'Detail {i}')
        
        plt.tight_layout()
        plt.show()

# 使用示例
denoiser = EEGDenoiser()
eeg_data = load_eeg_data('Z001.txt')
clean_data = denoiser.denoise(eeg_data)
denoiser.plot_decomposition(eeg_data)

8. 进阶技巧与优化建议

1. 自适应阈值选择：

   def adaptive_threshold(coeffs):
       thresholds = []
       for i in range(1, len(coeffs)):
           sigma = np.median(np.abs(coeffs[i])) / 0.6745
           thresholds.append(sigma * np.sqrt(2*np.log(len(coeffs[i]))))
       return thresholds
   

2. 多小波融合：结合不同小波基的优势

   def multi_wavelet_denoise(signal, wavelets=['db4', 'sym4'], level=4):
       results = []
       for w in wavelets:
           results.append(wavelet_denoise(signal, w, level))
       return np.mean(results, axis=0)
   

3. 实时处理优化：使用离散小波变换（DWT）的滑动窗口实现

在实际EEG分析中，我发现db4小波在大多数情况下表现稳定，但对某些特定频段的肌电干扰，结合sym4小波能获得更好的效果。另外，阈值的选择需要根据具体设备噪声特性进行调整，通常需要10-20组数据的测试来确定最佳参数。