脑电信号基于matlab小波变换DWT脑电信号ECG去噪含Matlab源码 1622期

Posted 2022-01-12 紫极神光

一、小波变换简介

大脑是由亿万个神经元组成的复杂系统，负责身体的各个功能的协调运作，通过大脑皮层上电极记录下的大脑细胞群的电位活动称为脑电信号。通过对脑电信号进行研究可以获得丰富的心理及生理的疾病信息，所以脑电信号的分析及去噪算法的研究无论是在临床的诊断和急病治疗上都是十分重要的。一般来说，脑电信号具有背景噪声强、信号微弱等特点，所以如何消除脑电数据的噪声，更好地获取对大脑有用的信息是当今研究的热门话题。

近年来小波变换广泛应用于信号处理和图像融合等领域。小波去噪的方法主要有模极大值法、空间滤波、阈值法，其中阈值法国内外学者广泛应用的方法[2]。阈值的选取是根据信号的长度设定，其设定的方法具有单一性，无法一次性对信号参数进行选定，由于选取过于零散，无法满足信号复杂等特点。

1 小波变换
小波变换是对基本小波进行伸缩和平移得到的，小波系数是原信号与小波基函数的相似系数。设函数Φ(t)为可积函数且满足t=L2(R），则傅里叶变换满足如下条件：

2 小波基的选择
小波基的选择对脑电去噪效果至关重要，一般对称性好的小波基不会产生相位畸变[8]。正则性好的小波基在进行脑电重构时脑电信号更加平滑，紧支撑的小波在脑电去噪时处理的速度更快。

本文选择的小波基为db Nv小波，可以满足以上的条件，db1（简称Haar）小波是小波基当中最简单的小波形式，局域性较差，因此本文选择db6小波，db6的特点是处理速度更快，MATLAB软件应用方便，其波形和脑电波形相似，对脑电去噪效果最好。db6小波曲线如图1所示。

图1 db6小波函数
3 小波去噪原理
小波去噪原理是将小波进行分解，再对小波进行多尺度变换，尽可能提取有用的脑电信号。然后再根据波恩脑电信号特征和噪声特点，选择合适的去噪模型，利用贝叶斯估计后的系数进行小波重构，从而得到去噪后的脑电信号。其基本原理如图2所示，在进行小波变换前，脑电信号的突变处峰值很高，突变处集中了噪声，很多平稳有用的脑电信号幅值很小，从而大大增加了去噪难度，所以本文以拉普拉斯为去噪模型，结合贝叶斯进行小波估计，已达到去噪的目的。

图2 脑电去噪原理

二、部分源代码

clear all;
%-----------------------------------
X = load ('118e00m_n.mat');
X = X.val;
%-----------------------------------
N = length(X);     %number of samples from original signal
z = nextpow2(N);
p = 2^z;
%update number samples if necessary by padding zeros at the end for a power
%of 2:
if(p > N)
    p = p-N;
    X(end + p)=0;
    N = length(X);
end
%-----------------------------------
%2^7=128, 2^8=256 2^9=512
%-----------------------------------
i_max = z;   %last level possible.  i = 0 is the first level
i=1;
%-----------------------------------
cD_cell = cell(1,1);  %cell unit of 1x1 dimension
cA_cell = cell(1,1);  %cell unit of 1x1 dimension
%-----------------------------------
%Decomposition
cA = X;
test=N;

while (i < i_max +1)
    [cA,cD]=dwt(cA,'coif4','mode','per');   %send current approximation and save cA and cD
    cD_celli,1= cD;
    cA_celli,1= cA;
    test = N/(2^i);
    i=i+1;
end

i=1;
while (i < i_max +1)
    cD_i = cD_celli,1;  %current vector of coefficient at level i

    l_d = length(cD_i);
    m_i = median(abs(cD_i));
    sig_i = m_i/0.6745;   %constant for threshold
    T = sig_i*sqrt(2.*log(l_d));
    k=1;
    while(k < l_d+1)
        if( abs(cD_i(k)) < T )
            cD_i(k)=0;
        end
            k=k+1;
    end
    cD_celli,1 = cD_i;
    i=i+1;
end
%-----------------------------------
%Reconstruction
i=i_max;
cA = cA_celli,1;
while (i > 0)
    cA = idwt(cA,cD_celli,1,'coif4','mode','per');
    i=i-1;
end
figure(2)
X = X(1:1024);
plot(X);
grid on;
title('Original Signal');
xlabel('Number of Samples');
ylabel('Amplitude in mV');

figure(3)
cA = cA(1:1024);
plot(cA);
grid on;
title('Denoised Signal');
xlabel('Number of Samples');
ylabel('Amplitude in mV');

%TODO: Refine algorithm for locating peaks..
% find R-peak every first peak

[~,locs_R]=findpeaks(X,'minpeakheight',50,'minpeakdistance',250);
X_inverted = -X;
[~,locs_S]=findpeaks(X_inverted,'minpeakheight',50,'minpeakdistance',250);

%--------------------------------------
[~,locs_R_D]=findpeaks(cA,'minpeakheight',50,'minpeakdistance',250);
cA_inverted = -cA;
[~,locs_S_D]=findpeaks(cA_inverted,'minpeakheight',50,'minpeakdistance',250);

%for detecting Q, from R-peak position find the minimum all the way back and save
%the spot
i= length(locs_R);
j=1;
locs_Q = zeros(length(i));

while(j< i+1)
a=locs_R(j)-35;
if(a<0)
    a=1;
end
temp = X_inverted(a:locs_R(j));   % size properly 

[~,index]=findpeaks(temp,'npeaks',1);
    b = length(temp)- index;
    locs_Q(j)= locs_R(j)-b;
    j=j+1;
end

三、运行结果

四、matlab版本及参考文献

1 matlab版本
2014a

2 参考文献
[1] 沈再阳.精通MATLAB信号处理[M].清华大学出版社，2015.
[2]高宝建,彭进业,王琳,潘建寿.信号与系统——使用MATLAB分析与实现[M].清华大学出版社，2020.
[3]王文光,魏少明,任欣.信号处理与系统分析的MATLAB实现[M].电子工业出版社，2018.
[4]徐洁.基于小波分析的脉搏波信号处理[J].电子设计工程. 2013,21(11)
[5]王宏旭,张晨洁,刘勇,郭滨.基于贝叶斯估计的小波脑电信号去噪算法研究[J].长春理工大学学报(自然科学版). 2021,44(02)