基于matlab的脑电信号处理
基于MATLAB的脑电信号处理研究

信息科学科技创新导报 Science and Technology Innovation Herald81脑科学的发展有着漫长的历史。
1924年,B e g e r首先记录并分析了脑电图,1933年,他的研究得到了英国著名生理学家E.D.Ad rian等人的肯定,从此脑电图得到迅速的发展,先后出现了多导联常规脑电图(自发脑电、诱发脑电)、睡眠脑电图、动态脑电图、视频脑电图等检测方法,研究的重点转移到了临床应用方面。
1 脑电信号的基础知识1.1 脑电信号与脑电图大脑皮层由150亿个各种形态的神经元及其复杂的纤维联系所组成。
神经元不仅在兴奋时伴有电位变化同时在抑制时或静止时也存在着电位现象。
我们把脑神经细胞的放电现象称为脑电,脑电活动是脑内大量神经元活动的综合表现。
脑电按频率划分主要有α波、β波、θ波、δ波等。
1.2 脑电信号的分类研究中的脑电信号可以分为两类:自发电位、诱发电位。
自发电位是一些特定的脑电成分,如,特殊皮层位置的特定脑电频率,它产生于正常脑功能过程中,与外部刺激无关;诱发电位是通过刺激产生的脑电信号,比如:由特定刺激诱发,通过潜伏期、幅度和位置识别的EEG波形。
2 相干平均法2.1 相关的基本概念(1)随机信号:在信号处理领域中,把信号分为确定性信号和随机性信号两大类型。
(2)噪声:通常人们处理信号是为了寻找某个特定的规律。
(3)非平稳过程:生物体具有适应外界环境的能力,当外界环境改变时,生物系统的参数会自动调整,信号的统计特性会随之改变,造成了生物医学信号的非平稳性质。
凡是不具备平稳随机过程概率属性的随机过程被称为非平稳随机过程。
(4)无偏和一致的估计:平稳信号统计特性的估计是生物医学信号处理的最基本的内容,其中评价估值的质量是一个很重要的工作。
2.2 相干平均法多个实测信号样本以时间基准点对齐,再将与同一时间对应的各样本数据求和平均,即可确定诱发响应的估值曲线,这种估计方法称为相干平均技术或同步平均技术。
如何使用Matlab进行神经电信号处理与分析

如何使用Matlab进行神经电信号处理与分析引言神经电信号是用于在神经系统中进行信息传递的电信号。
对于神经科学研究人员来说,了解和分析神经电信号对于揭示大脑的功能和机制至关重要。
Matlab作为一种功能强大的编程语言和开发环境,提供了丰富的工具和函数,可用于处理和分析神经电信号。
本文将介绍如何使用Matlab进行神经电信号处理和分析的基本方法和技巧。
1. 数据导入与预处理首先,我们需要将神经电信号数据导入Matlab环境中进行处理和分析。
常见的神经电信号数据包括脑电图(EEG)、脑磁图(MEG)和神经肌电图(EMG)。
Matlab提供了一系列函数和工具箱来读取不同格式的神经电信号数据文件,如EEGLAB和FieldTrip等。
在导入数据之后,我们还需要对数据进行预处理,以去除噪声和伪装信号。
常见的预处理步骤包括滤波、伪迹去除和去除运动伪迹等。
Matlab提供了诸多滤波函数和技术,如带通滤波、低通滤波和高通滤波等。
通过合理运用这些工具,我们可以对神经电信号数据进行有效的预处理。
2. 特征提取与分析在预处理之后,我们可以开始对神经电信号数据进行特征提取和分析。
特征提取是神经电信号处理中的重要步骤,可以从原始信号中提取出有用的信息。
在Matlab中,我们可以使用一系列函数和算法来提取信号的频域特征、时域特征和空域特征。
频域特征可以用于分析信号的频率成分和频率特征,如功率谱密度和频谱特性等。
Matlab提供了FFT和Periodogram等函数,可用于计算信号的频谱与谱密度。
时域特征包括信号的能量、振幅、时延和相位等。
可以使用Matlab的瞬时特征函数来计算信号的瞬时属性。
空域特征用于分析信号在空间上的分布和变化,如空间滤波和空间相关性等。
Matlab提供了一些函数和工具箱,如EEGLAB和FieldTrip,可用于实现这些分析。
3. 数据可视化与结果展示在特征提取和分析之后,我们通常需要将结果进行可视化展示,以便更好地理解和传达。
脑电gamma频段能量matlab

脑电gamma频段能量matlab1. 脑电信号及gamma频段人脑的神经元通过电信号的方式传递信息,这些电信号可以通过脑电图(Electroencephalography, EEG)来记录。
EEG信号通常被分类为不同的频段,例如delta、theta、alpha、beta和gamma频段,其中gamma频段通常指30-100Hz的频率范围。
近年来,研究表明gamma频段是大脑的一种重要信号,与人类的大脑功能密切相关,例如视觉注意力、知觉、记忆、思考等。
2. 脑电信号处理EEG信号是一种生理信号,通常具有复杂的时间关系和波形特征。
因此,对这种信号进行进一步分析需要采用一些特殊的方法。
Matlab是一种广泛使用的处理EEG信号的工具,其中一些常见的处理步骤包括:(1)预处理:通常包括去除各种干扰和噪声、滤波和对信号进行归一化等。
(2)信号分析:应用各种技术来识别、分类和研究EEG信号的不同频段和波形形状等。
(3)可视化:通过图形界面或者其他方法来呈现EEG信号,方便研究人员对其进行进一步分析和研究。
对于一个EEG信号,通常需要计算其各个频段的能量(power spectrum),以便研究者可以了解该信号在不同频段上的激活程度。
对于gamma频段,通常使用以下步骤来计算其能量:(1)首先,将原始EEG信号滤波,只保留目标频段(例如30-100Hz)内的信号。
(2)然后,将滤波后的信号平方,以获得其能量。
(3)最后,对于每一个时间点,将一段时间内的信号能量进行求平均,以得到该时间段内的平均gamma频段能量。
以下是在Matlab中完成这个任务的一些示例代码:% 假设我们已经有了原始EEG信号eegData和采样频率fs% 设定gamma频段下限和上限gammaLowLimit = 30;gammaHighLimit = 100;% 滤波gamma频段内的信号[b, a] = butter(4, [gammaLowLimit/(fs/2) gammaHighLimit/(fs/2)]);gammaFilteredData = filtfilt(b, a, eegData);% 计算信号能量gammaPower = gammaFilteredData.^2;在这个代码中,我们首先使用butter函数设计一个Butterworth低通滤波器,以保留30-100Hz的gamma频段内的信号。
matlab脑电信号处理

matlab脑电信号处理matlab脑电信号处理t=0.001:0.001:1;x=load('C:\Users\yxzhang\Desktop\rest_close.txt'); %读取文件y=load('C:\Users\yxzhang\Desktop\audio_close.txt');xx={}; %每个导联的数据存储yy={};n=1000; %数据数目sc=7; %小波包的分解尺度for i=1:1:8 %导联的数据分离xx{i}=x(:,i);yy{i}=y(:,i);endfor i=1:1:8%画出原始信号图像figuresubplot(2,2,1)plot(t,xx{i})axis([min(t) max(t) 1.1*floor(min(xx{i})) 1.1*ceil(max(xx{i}))])title('rest close 原始信号')ylabel('幅值')subplot(2,2,2)plot(t,yy{i})axis([min(t) max(t) 1.1*floor(min(yy{i})) 1.1*ceil(max(yy{i}))])title('audio close 原始信号')ylabel('幅值')%fft_原始信号的频谱分析xx1=fft(xx{i},n);pxx1=xx1.*conj(xx1)/n;yy1=fft(yy{i},n);pyy1=yy1.*conj(yy1)/n;%画出0-30hz内的功率谱图像n=60;subplot(2,2,3);plot(f,pxx1(1:n/2));ylabel('功率谱幅值(mv^2)');title('rest_close功率谱')subplot(2,2,4);plot(f,pyy1(1:n/2));ylabel('功率谱幅值(mv^2)');title('audio_close的功率谱')%分解信号选择分解尺度为7,同时重构信号wpt=wpdec(xx{i},sc,'db7','shannon'); %小波包分解信号xx80=wprcoef(wpt,[sc,0]); %重构信号xx81=wprcoef(wpt,[sc,1]);xx82=wprcoef(wpt,[sc,2]);xx83=wprcoef(wpt,[sc,3]);wpt=wpdec(yy{i},8,'db7','shannon');yy80=wprcoef(wpt,[sc,0]);yy81=wprcoef(wpt,[sc,1]);yy82=wprcoef(wpt,[sc,2]);yy83=wprcoef(wpt,[sc,3]);%画出重构信号figuresubplot(2,1,2);plot(yy80);title('audio close delta');ylabel('幅值');subplot(2,1,1);plot(xx80);title('rest close delta');ylabel('幅值');figuresubplot(2,1,1);plot(xx81);title('rest close theta');ylabel('幅值');subplot(2,1,2);plot(yy81);title('audio close theta');ylabel('幅值');subplot(2,1,1);plot(xx82);title('rest close alpha'); ylabel('幅值');subplot(2,1,2);plot(yy82);title('audio close alphta'); ylabel('幅值');figuresubplot(2,1,1);plot(xx83);title('rest close beta');ylabel('幅值')subplot(2,1,2);plot(yy83);title('audio close beta'); ylabel('幅值');n=1000;%fft_重构信号的频谱分析xx180=fft(xx80,n);pxx180=xx180.*conj(xx180)/n;xx181=fft(xx81,n);pxx181=xx181.*conj(xx181)/n;xx182=fft(xx82,n);pxx182=xx182.*conj(xx182)/n;xx183=fft(xx83,n);pxx183=xx183.*conj(xx183)/n;yy180=fft(yy80,n);pyy180=yy180.*conj(yy180)/n;yy181=fft(yy81,n);pyy181=yy181.*conj(yy181)/n;yy182=fft(yy82,n);pyy182=yy182.*conj(yy182)/n;yy183=fft(yy83,n);pyy183=yy183.*conj(yy183)/n;%画出重构信号的功率谱图n=60;f=1:n/2;figuresubplot(2,1,1);plot(f,pyy180(1:n/2));ylabel('功率谱幅值(mv^2)');title('rest close delta的功率谱'); subplot(2,1,2);plot(f,pyy180(1:n/2));ylabel('功率谱幅值(mv^2)'); title('audio close detta的功率谱'); figuresubplot(2,1,1);plot(f,pxx181(1:n/2));ylabel('功率谱幅值(mv^2)'); title('rest close theta的功率谱'); subplot(2,1,2);plot(f,pyy181(1:n/2));ylabel('功率谱幅值(mv^2)'); title('audio close theta的功率谱'); figuresubplot(2,1,1);plot(f,pxx182(1:n/2));ylabel('功率谱幅值(mv^2)'); title('rest close alpha的功率谱'); subplot(2,1,2);plot(f,pyy182(1:n/2));ylabel('功率谱幅值(mv^2)'); title('audio close alpha的功率谱'); figuresubplot(2,1,1);plot(f,pxx183(1:n/2));ylabel('功率谱幅值(mv^2)'); title('rest close beta的功率谱'); subplot(2,1,2);plot(f,pyy183(1:n/2));ylabel('功率谱幅值(mv^2)'); title('audio close beta的功率谱'); end。
基于Matlab的脑电波信号处理

做脑电波信号处理滴嘿嘿。
Matlab addictedCodes%FEATURE EXTRACTERfunction [features] = EEGfeaturetrainmod(filename,m)a = 4;b = 7;d = 12;e = 30;signals = 0;for index = 1:9; % read in the first ten EEG data because the files are numbered as ha11test01 rather than ha11test1.s = [filename '0' num2str(index) '.dat'];signal = tread_wfdb(s);if signals == 0;signals = signal;else signals = [signals signal];endendfor index = 10:1:m/2; % read in the rest of the EEG training datas = [filename num2str(index) '.dat'];signal = tread_wfdb(s);if signals == 0;signals = signal;else signals = [signals signal];endend%%%%% modification just for varying the training testing ratio ------for index = 25:1:25+m/2; % read in the rest of the EEG training data s = [filename num2str(index) '.dat'];signal = tread_wfdb(s);if signals == 0;signals = signal;else signals = [signals signal];endend%%%%%end of modification just for varying the training testing ratio-----for l = 1:m % exrating features (power of each kind of EEG wave forms) [Pxx,f]=pwelch(signals(:,l)-mean(signals(:,l)), [], [], [], 200); % relative power fdelta(l) = sum(Pxx(find(f<a)));ftheta(l) = sum(Pxx(find(f<b&f>a)));falpha(l) = sum(Pxx(find(f<d&f>b)));fbeta(l) = sum(Pxx(find(f<e&f>d)));fgama(l)= sum(Pxx(find(f>e))); % gama wave included for additional workendfeatures = [fdelta; ftheta; falpha; fbeta a; fgama];features = features';end%CLASSIFIER%(Has three similar classification modifation: EEGclassification, EEGclassificationmod and EEGclassificationmod1 saved and used in the running file for additional works)function [class, err, classall, errall]= EEGclassification(trainfilename, m, testfilename,n, p,q)% p - waveform 1, q - wave form two o –wave form three% 1 - delta 2 - theta 3 - alpha 4 –beta 5 - Gamma[featurestrain] = EEGfeature(trainfilename, m);% modification to EEGfeaturemod function for varying testing training ratio[featurestest] = EEGfeature(testfilename,n);training = [featurestrain(:,p) featurestrain(:,q)];% modify how many features to extract heresample = [featurestest(:,p) featurestest(:,q)];group = [ones(m/2,1);2*ones(m/2,1)];traininga = featurestrain;samplea = featurestest;[class, err, POSTERIOR, logp, coeff]= classify(sample, training, group, 'quadratic'); %'mahalanobis','quadratic','linear'as default[classall, errall]= classify(samplea, traininga, group, 'quadratic');display(class);display(err);% running file%------------------ using 2 features out of 4 comparison -----------class = 0; err = 0;p = 1;for q = 2:1:4[clas, er]= EEGclassification('ha11train',50,'ha11test', 10, p,q);if class == 0; % appending newly generated classificaton and errorclass = clas;else class = [class clas];endif err ==0;err = er;else err = [err er];endendp = 2;for q = 3:4[clas, er]= EEGclassification('ha11train',50,'ha11test',10, p,q);class = [class clas];err = [err er]; % appending newly generated classificaton and error endp=3;q=4;[clas, er,classall, errall]= EEGclassification('ha11train',50,'ha11test', 10, p,q);class = [class clas classall];err = [err er errall];results = [class;err]; % displaying the results in a tabledisplay (results);%------------------ using 2 features out of 5...。
Matlab技术在人体生理信号处理中的应用案例解析

Matlab技术在人体生理信号处理中的应用案例解析引言:在当今数字化的时代,人体生理信号的处理越来越重要。
准确地捕捉和分析这些信号有助于我们深入了解人体的生理状况和健康状况。
而Matlab作为一款强大的工具,被广泛应用于人体生理信号的处理。
本文将通过几个应用案例,来展示Matlab技术在人体生理信号处理中的重要作用。
一、脑电图(EEG)处理脑电图是记录大脑电活动的一种生理信号。
它反映了大脑在不同状态下的电活动,被广泛用于研究睡眠、意识状态和神经疾病等。
在脑电图处理中,Matlab提供了丰富的工具和函数来处理和分析原始数据,例如滤波、频域分析和时间域分析等。
通过Matlab的工具箱,可以方便地对脑电图进行去除噪声、频谱分析、时频分析等,进而提取出感兴趣的特征。
这些特征可以用于研究不同的脑电图模式,并为医生提供诊断依据。
二、心电图(ECG)处理心电图是记录心脏电活动的一种生理信号。
它被广泛用于心脏疾病的诊断和监测。
Matlab提供了一套丰富的工具箱,例如信号处理工具箱和心电图工具箱,用于对心电图进行处理和分析。
通过Matlab的滤波技术和傅里叶变换等算法,可以去除心电图中的噪声和干扰,同时还可以提取心电图中的特征,例如心率、ST段变异和QRS波形等。
这些特征对于心脏疾病的早期诊断和监测非常重要。
三、肌电图(EMG)处理肌电图是记录肌肉电活动的一种生理信号。
它可以反映肌肉的收缩和放松状况,被广泛用于研究肌肉活动和肌肉疾病。
在肌电图处理中,Matlab提供了一系列用于去除噪声和分析肌电图的函数和工具箱。
通过Matlab的滤波算法和傅里叶变换等技术,可以准确地提取肌电图中的肌肉活动特征,并用于研究肌肉疾病和运动控制等领域。
四、血氧饱和度(SpO2)处理血氧饱和度是衡量氧气输送至组织的重要指标,通常通过红外线光源和传感器来测量。
Matlab提供了一些工具箱用于对血氧饱和度进行监测和处理。
通过Matlab的数据处理和曲线拟合技术,可以对原始的SpO2数据进行滤波和校正,提高数据的准确性和可信度。
Matlab中的脑电图信号处理与频谱分析方法

Matlab中的脑电图信号处理与频谱分析方法一、引言脑电图(Electroencephalogram,简称EEG)是记录大脑电生理活动的一种非侵入性方法。
在临床和研究中,脑电图被广泛应用于诊断神经系统疾病、研究认知过程等领域。
而Matlab作为一种功能强大的科学计算软件,提供了丰富的工具箱和函数,可以有效地处理和分析脑电图数据。
本文将介绍Matlab中常用的脑电图信号处理方法和频谱分析技术。
二、脑电图信号处理方法1. 清除噪声在进行脑电图信号分析之前,首先需要对原始信号进行预处理,以去除噪声和伪迹。
Matlab提供了多种滤波器函数,如低通滤波器和带通滤波器,可以有效地去除不需要的高频噪声和低频噪声。
2. 分段处理脑电图信号常常是一个连续的时间序列,在某些情况下,可以将信号分成较短的时间段。
这样做有助于分析信号在不同时间段的特性。
Matlab中可以使用窗函数对信号进行分段处理,并通过遍历每个窗口进行连续的分析。
3. 时域分析时域分析是对信号在时间上的变化进行定量描述的方法。
常用的时域分析方法包括计算信号的平均值、方差、峰值和时域波形图等。
在Matlab中,可以使用相应的函数和工具箱进行时域特征提取和可视化,从而实现对脑电图信号的时域分析。
4. 频域分析频域分析是对信号在频率上的变化进行研究和描述的方法。
脑电图信号通常包含不同频率的成分,因此频域分析对于理解信号的特征和性质非常重要。
Matlab提供了多种频谱分析方法,如快速傅里叶变换(FFT)、小波变换(Wavelet Transform)和自相关函数等。
这些函数可以帮助我们从频域的角度来研究脑电图信号,并提取出频率成分的信息。
三、频谱分析方法1. 快速傅里叶变换(FFT)FFT是一种常见的频谱分析方法,可以将信号从时域转换为频域。
通过计算信号的幅度谱和相位谱,我们可以获得信号在不同频率下的能量分布。
Matlab 中的fft函数可以高效地计算快速傅里叶变换,并绘制出脑电图信号的频谱图。
基于matlab的脑电信号处理

基于matlab的脑电信号处理————————————————————————————————作者:————————————————————————————————日期:南京航空航天大学基于Matlab的脑电信号处理姓名陆想想专业领域生物医学工程课程名称数字信号处理二О一三年四月摘要:脑电信号属于非平稳随机信号,且易受到各种噪声干扰。
本文基于Matlab仿真系统,主要研究了小波变换在脑电信号处理方面的应用,包括小波变换自动阈值去噪处理、强制去噪处理,以α波为例,提取小波分解得到的各层频率段的信号,并做了一定的分析和评价。
关键词:脑电信号;小波变换;去噪重构;频谱分析0 引言脑电信号EEG(Electroencephalograph)是人体一种基本生理信号,蕴涵着丰富的生理、心理及病理信息,脑电信号的分析及处理无论是在临床上对一些脑疾病的诊断和治疗,还是在脑认知科学研究领域都是十分重要的。
由于脑电信号的非平稳性且极易受到各种噪声干扰,特别是工频干扰。
因此消除原始脑电数据中的噪声,更好地获取反映大脑活动和状态的有用信息是进行脑电分析的一个重要前提。
本文的研究目的是利用脑电采集仪器获得的脑电信号,利用Fourier变换、小波变换等方法对脑电信号进行分析处理,以提取脑电信号α波的“梭形”节律,并对脑电信号进行功率谱分析和去噪重构。
1 实验原理和方法1.1实验原理1.1.1脑电信号根据频率和振幅的不同,可以将脑电波分为4种基本类型[1],即δ波、θ波、α波、β波。
4种波形的起源和功能也不相同,如图1所示。
图1 脑电图的四种基本波形α波的频率为8~13Hz,振幅为为20~100µV,它是节律性脑电波中最明显的波,整个皮层均可产生α波。
正常成人在清醒、安静、闭目时,波幅呈现有小变大,再由大变小,如此反复进行,形成所谓α节律的“梭形”。
每一“梭形”持续时间约为1~2s。
当被试者睁眼、警觉、思考问题或接受其他刺激时,α波立即消失而代之以快波,这种现象称之为“α波阻断”。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
基于matlab的脑电信号处理————————————————————————————————作者:————————————————————————————————日期:南京航空航天大学基于Matlab的脑电信号处理姓名陆想想专业领域生物医学工程课程名称数字信号处理二О一三年四月摘要:脑电信号属于非平稳随机信号,且易受到各种噪声干扰。
本文基于Matlab仿真系统,主要研究了小波变换在脑电信号处理方面的应用,包括小波变换自动阈值去噪处理、强制去噪处理,以α波为例,提取小波分解得到的各层频率段的信号,并做了一定的分析和评价。
关键词:脑电信号;小波变换;去噪重构;频谱分析0 引言脑电信号EEG(Electroencephalograph)是人体一种基本生理信号,蕴涵着丰富的生理、心理及病理信息,脑电信号的分析及处理无论是在临床上对一些脑疾病的诊断和治疗,还是在脑认知科学研究领域都是十分重要的。
由于脑电信号的非平稳性且极易受到各种噪声干扰,特别是工频干扰。
因此消除原始脑电数据中的噪声,更好地获取反映大脑活动和状态的有用信息是进行脑电分析的一个重要前提。
本文的研究目的是利用脑电采集仪器获得的脑电信号,利用Fourier变换、小波变换等方法对脑电信号进行分析处理,以提取脑电信号α波的“梭形”节律,并对脑电信号进行功率谱分析和去噪重构。
1 实验原理和方法1.1实验原理1.1.1脑电信号根据频率和振幅的不同,可以将脑电波分为4种基本类型[1],即δ波、θ波、α波、β波。
4种波形的起源和功能也不相同,如图1所示。
图1 脑电图的四种基本波形α波的频率为8~13Hz,振幅为为20~100µV,它是节律性脑电波中最明显的波,整个皮层均可产生α波。
正常成人在清醒、安静、闭目时,波幅呈现有小变大,再由大变小,如此反复进行,形成所谓α节律的“梭形”。
每一“梭形”持续时间约为1~2s。
当被试者睁眼、警觉、思考问题或接受其他刺激时,α波立即消失而代之以快波,这种现象称之为“α波阻断”。
一般认为,α波是大脑皮质处于清醒安静状态时电活动的主要表现。
β波的频率是18~30Hz,振幅为5~20µV,是一种快波。
β波的出现以为着大脑比较兴奋。
θ波的频率是4~7Hz,振幅为10~50µV,它是在困倦时,神经系统处于抑制状态时所记录的波形。
δ波在睡眠、深度麻醉、缺氧或大脑有器质性病变时出现,频率是1~3.5Hz,振幅为20~200µV。
1.1.2小波变换小波变换的概念是由从事石油信号处理的法国工程师J.Morlet在1974年首先提出的。
与Fourier变换相比,小波变换是时间(空间)频率的局部化分析,它通过伸缩平移运算,对信号逐步进行多尺度细化,最终达到高频处时间细分,低频处频率细分,能自动适应时频信号分析的要求,从而可以聚焦到信号的任意细节,解决了Fourier变换的困难问题。
在噪声中如何准确地检测到信号一直是信号处理领域所关心的内容。
小波变换是一种信号的时间一尺度分析方法,由于具有多分辨率分析的特点,良好的时频局部化特性,能够对各种时变信号进行有效的分解,从而较好地将信号与噪声加以分离,获得满意的去噪效果。
小波分析去噪处理的方法一般有三种:默认阈值去噪处理。
该方法是利用ddencmp函数生成信号的默认阈值,然后利用wdencmp函数进行去噪处理。
给定阈值进行去噪处理。
在实际的去噪过程中,阈值往往可以通过经验公式获得,并且这种阈值比默认阈值可信度高。
在进行量化处理时可利用wthresh函数。
强制去噪处理。
该方法是将小波分解结构中的高频或者低频系数设置为0,即滤掉所有高频部分或低频部分。
这种方法比较简单,且去噪得到的信号比较平滑,但是容易丢失信号中的有用成分[2]。
本文采用了两种去噪方法,并分析比较了他们的去噪效果。
1.2实验方法与步骤1.2.1脑电信号的读取本文使用的脑电数据是使用南航生物医学光子学实验室的脑电采集系统采集获得的,原始数据格式为.eeg。
为了方便在Matlab环境下对数据进行分析,将文件转换为.txt格式。
脑电采集使用的是16通道,采样频率为256Hz,文件中存储的数据的形式为数据点数×通道数。
实验中选取了第14通道的前8000个数据点作为样本进行分析。
由于采样时间是256Hz,所以这段信号的持续时间大约是32秒。
1.2.2信号的频域和功率谱分析为了研究脑电信号中不同频率信号的能量分布以及变化情况,首先对样本信号进行Fourier变换,得到频域图。
然后进一步对信号做功率谱分析,得到功率谱图,从功率谱图中,可以直观的观察到不同频率信号的能量分布情况。
由于脑电数据是在被采集者安静清醒的状态下采集得到的,理论上α波应该占主导地位。
1.2.3信号的小波变换及重构基于小波变换降噪处理的方法通常有3个步骤:首先是将信号进行n层小波变换,得到小波系数;然后在小波变换域上利用信号与噪声的不同特性,对小波变换进行阈值化处理,把噪声从信号中区分开来(主要是对高频系数进行阈值化处理;最后是利用重构算法重构信号。
小波变换去噪的效果主要取决于对含噪信号的噪声估计方法以及所采用的小波函数[4]。
本文使用其中两种去噪方法。
第一种是默认阈值去噪,先对样本信号进行8层小波分解,使用的小波函数是db8,然后利用ddencmp函数生成信号的默认阈值,最后利用wdencmp 函数进行去噪处理,得到去噪的信号。
另一种是强制去噪,同样地,先对样本信号进行8层小波分解,采用的小波函数也是db8。
然后提取各层小波系数,再将脑电信号频率以外的几层小波系数置零,得到重构的各层小波系数。
最后由重构的小波系数得到重构的信号。
该信号中除去了某些频率的信号,起到了去噪的效果。
2仿真结果及分析2.1原始脑电数据的读取和显示采集到的脑电信号文件为data.txt,调用eeg_load.m文件,即可绘制出脑电样本信号图,如下图2所示。
图2 脑电样本信号图2.2脑电信号频谱图及功率谱图的绘制首先调用eeg_fft.m文件,原理是对样本信号进行傅立叶变换[3],即可获得样本信号的频谱图,如下图3所示。
图3信号频谱图从频谱图中可以看出,低频信号和11Hz左右的信号特别强,25Hz以上的信号几乎为零。
由于α波的频率为8~13Hz,由此可以知道,该信号中α波比较多。
在脑电采集过程中,被采集者没有处于深度睡眠状态,因此接近0Hz的低频信号可以确定为噪声。
15~25Hz频段的信号很微弱,因此可以以判断出,信号中几乎没有β波。
谱估计是数字信号处理的主要内容之一,主要研究信号在频域中的各种特征,目的是根据有限数据在频域内提取被淹没在噪声中的有用信号。
调用eeg_psd,m文件,可以绘制信号的功率谱图,如下图4所示。
图4 信号功率谱图原始的功率谱图是关于f=128Hz对称的,为了便于分析,截取了0~65Hz的一段。
信号的功率谱显示,在10Hz作用频率处能量出现峰值,表明在平静状态下采集的脑电信号中α波能量最大,符合生理学的研究结论。
2.3脑电信号小波分解各层重构波形实验选取db8小波对前述采集的样本信号进行分解,调用wavelet_dec.m文件,分解得到的各层信号如下图5所示。
图5 脑电信号经小波分解后的各层分量由于采样频率是256Hz,所以对于自带的频率范围是:]2,2[],2,2[],2,2[],2,2[],2,2[],2,2[],2,2[],2,0[22334455667788s s s s s s s s s s s s s s s f f f f f f f f f f f f f f f相应的各子带的频率成分如下表所示:子带信号 频率范围 主要信号成分 a8 0Hz~1Hz 低频干扰 d8 1Hz~2Hz δ波 d7 2Hz~4H δ波 d6 4Hz~8Hz θ波 d5 8Hz~16Hz α波 d4 16Hz~32Hz β波 d3 32Hz~64Hz 高频噪声 d2 64Hz~128Hz 高频噪声 d1128Hz~256Hz高频噪声表1 小波分解后脑电信号子带的频率范围一般情况下,a8子带内是低频干扰,d3、d2、d1子带内是高频噪声,d8、d7、d6、d5、d4子带内是脑电信号,但也可能混有一定的噪声,需要根据实际情况来分析。
2.4 脑电信号节律提取α波是节律性脑电波中最明显的波。
由以上小波分解得到的各层分量分析可知,α波主要集中在d5子带内。
所以选取了d5子带内第5500~8000之间的2500个数据点,将信号图绘制出来,如下图6所示。
图6 提取出来的α波节律从上图中可以明显看出α波波幅由小到大,再由大到小作规律性变化,呈棱状图形。
2.5 小波分解去噪和重构波形2.5.1小波变换默认阈值去噪处理调用noise_reduce.m文件,可以实现脑电信号的小波变换默认阈值去噪处理,原始数据及去噪处理结果对比如下图7所示。
图7 原始信号与小波默认阈值去噪结果图的对比从原始脑电信号与去噪处理后的效果来看,经去噪处理后的信号高频信号有所减少。
2.5.2小波变换强制去噪处理调用wavelet_rec.m文件,绘制小波变换强制去噪处理之后的信号如下图8所示。
图8 原始信号与小波强制去噪结果图的对比该图表示的是除去了低频干扰和高频噪声之后的结果,从与原始信号的对比中可以看出,高频噪声很明显被消除了。
但由于考虑到实际情况,16~32Hz子带内极少有β波,大部分为噪声,所以把这个频带内的信号也全部清除了。
3 讨论基于FFT变换的功率谱估计适用于平稳时间信号分析,计算结果只能反映信号段总体的平均功率分布情况, 不包含信号的任何时域变化信息, 并且谱估计的频率分辨率与所采用的信号长度成正比,即受数据点数目的影响。
基于小波变换的去噪方法,对非平稳信号去噪,要比传统的滤波去噪声得到的效果好,主要是由于传统的滤波器都具有低通性,对需要分析在每个时刻含有不同频率成分的非平稳信号来说,是很难对它进行匹配分析。
而小波变换具有多分辨率,且在时频域都具有局部性,因此很适合用来分析非平稳信号。
在用小波分析来进行去噪的关键在于阈值的选取,如果阈值选取的太高,会使信号失去太多细节,使信号失真,如果阈值选得太低,又不能达到去噪的目的。
4 结论本文利用实测的原始脑电信号, 对脑电信号的处理方法与结果进行了一定的分析和评价, 以期为脑电信号处理及特征提取提供一定的理论参考和分析依据。
目前人们也尝试用非线性处理方法、神经网络的方法、时频结合等等现代的方法来处理脑电信号, 相信这些方法会为脑认知以及医学的发展作出贡献[5]。
脑电信号属于非平稳随机信号,小波分析的方法可以直接对信号的某些频率分量进行观察或者提取出有用的特征信号,为脑电信号的测量与分析提供了非常好的前景。