小波自适应去噪在脑电信号处理MATLAB仿真实现
1. 脑电信号特点与小波去噪原理
1.1 脑电信号特性
- 频率范围：0.5-100 Hz，主要能量集中在0.5-30 Hz
- 信号幅度：10-100 μV，非常微弱
- 噪声来源：
- 工频干扰：50/60 Hz电源干扰
- 肌电干扰：肌肉活动产生的高频噪声
- 眼电干扰：眨眼和眼球运动
- 基线漂移：低频成分
- 测量噪声：仪器本身噪声
1.2 小波去噪优势
- 时频局部化：能同时分析信号的时域和频域特征
- 多分辨率分析：适合分析非平稳的脑电信号
- 自适应特性：不同尺度采用不同阈值策略
2. 小波自适应去噪算法原理
2.1 基本流程
原始脑电信号 → 小波分解 → 系数阈值处理 → 小波重构 → 去噪后信号
2.2 关键算法步骤
步骤1：小波分解
选择合适的小波基函数和分解层数：
\([cA_n, cD_n, cD_{n-1}, ..., cD_1] = \text{DWT}(EEG)\)
其中：
- \(cA_n\)：第n层近似系数（低频）
- \(cD_k\)：第k层细节系数（高频）
步骤2：阈值选择策略
1. 通用阈值（VisuShrink）
\(\lambda = \sigma \sqrt{2\ln(N)}\)
其中(\sigma)为噪声标准差，(N)为信号长度。
2. 自适应阈值（SURE Shrink）
基于Stein无偏风险估计，最小化均方误差。
3. 分层阈值
不同分解层采用不同阈值：
\(\lambda_j = \frac{\sigma_j \sqrt{2\ln(N)}}{\sqrt{j}}\)
步骤3：阈值函数
硬阈值函数：
软阈值函数：
3. MATLAB仿真实现
3.1 完整的自适应去噪算法
%% 小波自适应去噪脑电信号仿真
clear; close all; clc;
%% 参数设置
fs = 256; % 采样频率(Hz)
t = 0:1/fs:10-1/fs; % 时间向量
N = length(t); % 信号长度
%% 生成模拟脑电信号（包含多种节律）
% 脑电节律成分
delta = 0.5 * sin(2*pi*2*t); % Delta波 (1-4 Hz)
theta = 0.8 * sin(2*pi*6*t); % Theta波 (4-8 Hz)
alpha = 1.2 * sin(2*pi*10*t); % Alpha波 (8-13 Hz)
beta = 0.6 * sin(2*pi*20*t); % Beta波 (13-30 Hz)
% 组合成纯净脑电信号
eeg_clean = delta + theta + alpha + beta;
%% 添加各种噪声
% 1. 高斯白噪声（仪器噪声）
noise_white = 0.3 * randn(1, N);
% 2. 工频干扰 (50Hz + 谐波)
noise_50hz = 0.5 * sin(2*pi*50*t) + 0.2 * sin(2*pi*100*t);
% 3. 肌电干扰 (突发性高频噪声)
emg_noise = zeros(1, N);
emg_positions = [1000, 1500, 2000, 3500, 4000]; % 突发噪声位置
for pos = emg_positions
duration = 50; % 持续50个采样点
idx = pos:min(pos+duration-1, N);
emg_noise(idx) = 1.5 * randn(1, length(idx));
end
% 4. 眼电干扰 (低频大幅值)
blink_noise = zeros(1, N);
blink_positions = [500, 1800, 3000];
for pos = blink_positions
duration = 100;
idx = pos:min(pos+duration-1, N);
blink_noise(idx) = 3 * (1 - cos(2*pi*(0:length(idx)-1)/duration));
end
% 含噪脑电信号
eeg_noisy = eeg_clean + noise_white + noise_50hz + emg_noise + blink_noise;
%% 小波自适应去噪函数
function eeg_denoised = wavelet_denoise_adaptive(eeg_signal, fs, wavelet_name, level)
% 小波自适应去噪
% 输入：eeg_signal - 含噪脑电信号
% fs - 采样频率
% wavelet_name - 小波名称
% level - 分解层数
% 输出：eeg_denoised - 去噪后信号
N = length(eeg_signal);
% 小波分解
[C, L] = wavedec(eeg_signal, level, wavelet_name);
% 提取各层系数
approx_coef = appcoef(C, L, wavelet_name, level);
detail_coefs = cell(level, 1);
for i = 1:level
detail_coefs{i} = detcoef(C, L, i);
end
% 自适应阈值处理
C_denoised = C; % 初始化去噪后系数
% 对细节系数进行分层阈值处理
start_idx = L(1) + 1;
for i = 1:level
coef_length = L(i+1);
current_coef = C(start_idx:start_idx+coef_length-1);
% 自适应阈值选择
if i <= 2 % 高频层使用更严格的阈值
threshold = thselect(current_coef, 'rigrsure'); % SURE阈值
else % 低频层使用较宽松的阈值
threshold = thselect(current_coef, 'heursure'); % 启发式阈值
end
% 软阈值处理（更好保留信号特征）
denoised_coef = wthresh(current_coef, 's', threshold);
% 更新系数
C_denoised(start_idx:start_idx+coef_length-1) = denoised_coef;
start_idx = start_idx + coef_length;
end
% 小波重构
eeg_denoised = waverec(C_denoised, L, wavelet_name);
% 确保长度一致
if length(eeg_denoised) > N
eeg_denoised = eeg_denoised(1:N);
elseif length(eeg_denoised) < N
eeg_denoised(end+1:N) = 0;
end
end
%% 执行小波去噪
wavelet_name = 'db4'; % Daubechies 4小波
level = 5; % 分解层数
eeg_denoised = wavelet_denoise_adaptive(eeg_noisy, fs, wavelet_name, level);
%% 性能评估
% 信噪比计算
function snr_val = calculate_snr(clean_signal, noisy_signal)
signal_power = mean(clean_signal.^2);
noise_power = mean((noisy_signal - clean_signal).^2);
snr_val = 10 * log10(signal_power / noise_power);
end
% 均方根误差
function rmse_val = calculate_rmse(signal1, signal2)
rmse_val = sqrt(mean((signal1 - signal2).^2));
end
% 计算评估指标
snr_input = calculate_snr(eeg_clean, eeg_noisy);
snr_output = calculate_snr(eeg_clean, eeg_denoised);
rmse_before = calculate_rmse(eeg_clean, eeg_noisy);
rmse_after = calculate_rmse(eeg_clean, eeg_denoised);
%% 时频分析函数
function [tfr, freq, time] = compute_stft(signal, fs, window_len, overlap)
% 短时傅里叶变换
window = hamming(window_len);
noverlap = round(overlap * window_len);
nfft = max(256, 2^nextpow2(window_len));
[S, F, T] = spectrogram(signal, window, noverlap, nfft, fs);
tfr = 20*log10(abs(S) + eps);
freq = F;
time = T;
end
%% 结果可视化
figure('Position', [100, 100, 1400, 1000]);
% 1. 时域信号对比
subplot(3,2,1);
plot(t, eeg_clean, 'b', 'LineWidth', 1.5); hold on;
plot(t, eeg_noisy, 'r', 'LineWidth', 0.5, 'Alpha', 0.7);
xlabel('时间 (s)'); ylabel('幅度 (\muV)');
title('原始脑电信号 vs 含噪信号');
legend('纯净信号', '含噪信号', 'Location', 'best');
grid on;
subplot(3,2,2);
plot(t, eeg_clean, 'b', 'LineWidth', 1.5); hold on;
plot(t, eeg_denoised, 'g', 'LineWidth', 1.2);
xlabel('时间 (s)'); ylabel('幅度 (\muV)');
title('原始脑电信号 vs 去噪后信号');
legend('纯净信号', '去噪信号', 'Location', 'best');
grid on;
% 2. 时频分析
window_len = 128; overlap = 0.75;
% 纯净信号时频图
subplot(3,2,3);
[tfr_clean, freq_clean, time_clean] = compute_stft(eeg_clean, fs, window_len, overlap);
imagesc(time_clean, freq_clean, tfr_clean);
axis xy; colorbar; clim([-20 40]);
xlabel('时间 (s)'); ylabel('频率 (Hz)');
title('纯净信号时频图');
% 含噪信号时频图
subplot(3,2,4);
[tfr_noisy, freq_noisy, time_noisy] = compute_stft(eeg_noisy, fs, window_len, overlap);
imagesc(time_noisy, freq_noisy, tfr_noisy);
axis xy; colorbar; clim([-20 40]);
xlabel('时间 (s)'); ylabel('频率 (Hz)');
title('含噪信号时频图');
% 去噪信号时频图
subplot(3,2,5);
[tfr_denoised, freq_denoised, time_denoised] = compute_stft(eeg_denoised, fs, window_len, overlap);
imagesc(time_denoised, freq_denoised, tfr_denoised);
axis xy; colorbar; clim([-20 40]);
xlabel('时间 (s)'); ylabel('频率 (Hz)');
title('去噪后信号时频图');
% 3. 性能指标显示
subplot(3,2,6);
text(0.1, 0.9, sprintf('输入信噪比: %.2f dB', snr_input), 'FontSize', 12);
text(0.1, 0.7, sprintf('输出信噪比: %.2f dB', snr_output), 'FontSize', 12);
text(0.1, 0.5, sprintf('SNR改善: %.2f dB', snr_output - snr_input), 'FontSize', 12);
text(0.1, 0.3, sprintf('RMSE (去噪前): %.4f', rmse_before), 'FontSize', 12);
text(0.1, 0.1, sprintf('RMSE (去噪后): %.4f', rmse_after), 'FontSize', 12);
axis off;
title('去噪性能指标');
%% 不同小波基函数比较
wavelets = {'db4', 'sym4', 'coif3', 'bior3.3'};
snr_improvements = zeros(1, length(wavelets));
rmse_results = zeros(1, length(wavelets));
fprintf('=== 不同小波基函数性能比较 ===\n');
for i = 1:length(wavelets)
eeg_temp = wavelet_denoise_adaptive(eeg_noisy, fs, wavelets{i}, level);
snr_temp = calculate_snr(eeg_clean, eeg_temp);
rmse_temp = calculate_rmse(eeg_clean, eeg_temp);
snr_improvements(i) = snr_temp - snr_input;
rmse_results(i) = rmse_temp;
fprintf('小波基: %s, SNR改善: %.2f dB, RMSE: %.4f\n', ...
wavelets{i}, snr_improvements(i), rmse_temp);
end
%% 绘制不同小波性能比较图
figure('Position', [200, 200, 1000, 400]);
subplot(1,2,1);
bar(snr_improvements, 'FaceColor', [0.2 0.6 0.8]);
set(gca, 'XTickLabel', wavelets);
ylabel('SNR改善 (dB)');
title('不同小波基函数的SNR改善');
grid on;
subplot(1,2,2);
bar(rmse_results, 'FaceColor', [0.8 0.4 0.2]);
set(gca, 'XTickLabel', wavelets);
ylabel('RMSE');
title('不同小波基函数的RMSE');
grid on;
3.2 高级自适应算法扩展
%% 基于小波包的自适应去噪（更精细的频率划分）
function eeg_denoised_wp = wavelet_packet_denoise(eeg_signal, wavelet_name, level)
% 小波包去噪 - 提供更精细的频带划分
% 小波包分解
wp = wpdec(eeg_signal, level, wavelet_name);
% 计算最佳基（使用熵准则）
entropy_type = 'shannon';
wp_best = bestbas(wp, entropy_type);
% 节点阈值处理
for i = 0:(2^level - 1)
node_coef = wpcoef(wp_best, [level, i]);
if ~isempty(node_coef)
% 自适应阈值
threshold = thselect(node_coef, 'rigrsure');
denoised_coef = wthresh(node_coef, 's', threshold);
% 更新系数
wp_best = write(wp_best, 'cfs', [level, i], denoised_coef);
end
end
% 小波包重构
eeg_denoised_wp = wprec(wp_best);
end
%% 执行小波包去噪
eeg_wp_denoised = wavelet_packet_denoise(eeg_noisy, 'db4', 4);
% 比较性能
snr_wp = calculate_snr(eeg_clean, eeg_wp_denoised);
fprintf('\n小波包去噪性能: SNR = %.2f dB\n', snr_wp);
4. 实际脑电信号处理应用
4.1 真实脑电数据处理建议
%% 实际脑电数据处理流程
function processed_eeg = process_real_eeg(eeg_data, fs, channels)
% 实际脑电信号处理流程
% 输入：eeg_data - 多通道脑电数据
% fs - 采样率
% channels - 要处理的通道索引
[num_channels, num_samples] = size(eeg_data);
processed_eeg = zeros(length(channels), num_samples);
for i = 1:length(channels)
ch_idx = channels(i);
raw_signal = eeg_data(ch_idx, :);
% 1. 预处理：去除基线漂移
baseline = medfilt1(raw_signal, fs); % 中值滤波估计基线
signal_detrend = raw_signal - baseline;
% 2. 工频陷波 (50Hz)
wo = 50/(fs/2); bw = wo/35;
[b,a] = iirnotch(wo, bw);
signal_notched = filtfilt(b, a, signal_detrend);
% 3. 小波自适应去噪
signal_denoised = wavelet_denoise_adaptive(signal_notched, fs, 'db4', 5);
processed_eeg(i, :) = signal_denoised;
end
end
参考代码 小波自适应去噪脑电信号 www.youwenfan.com/contentcnl/80597.html
5. 算法性能分析
5.1 优势
- 自适应性强：根据不同频带特性自动调整阈值
- 保留特征：软阈值更好地保留脑电信号特征波形
- 多分辨率：同时处理高频噪声和低频干扰
- 计算效率：相比传统滤波方法有更好时频局部化
5.2 参数选择建议
- 小波基：推荐
db4、sym4，平衡光滑性和紧支撑 - 分解层数：4-6层，根据采样率调整
- 阈值策略：高频层使用严格阈值，低频层使用宽松阈值
- 阈值函数：软阈值更适合脑电信号
5.3 临床应用价值
- 癫痫检测：有效去除噪声，突出棘波、尖波特征
- 睡眠分期：清晰分离不同睡眠阶段的脑电节律
- 脑机接口：提高信号质量，提升分类准确率
- 认知研究：更好地分析事件相关电位
这个完整的小波自适应去噪框架为脑电信号处理提供了强大的工具，能够有效应对脑电信号中的各种噪声干扰，同时很好地保留有用的生理信息特征。
浙公网安备 33010602011771号