脑电数据处理分析——edf转mat及fft

最新推荐文章于 2024-10-24 16:41:31 发布

hi_linda

最新推荐文章于 2024-10-24 16:41:31 发布

阅读量7.3k

点赞数 9

分类专栏：脑电文章标签：脑电fft

本文链接：https://blog.csdn.net/qq_23100417/article/details/108142536

版权

脑电专栏收录该内容

15 篇文章

订阅专栏

声明：参考来源：https://zhuanlan.zhihu.com/p/34410111

生理信号数据常用的格式为edf，而matlab常用的是mat,txt。

如果已经在数据库里下载了edf格式的数据。不妨先了解下EDF格式具体内容：

HEADER RECORD

8 ascii : 数据格式的版本
80 ascii : 被试ID
80 ascii : 数据记录编号

8 ascii : 开始记录的日期 (dd.mm.yy)
8 ascii : 开始记录的时间 (hh.mm.ss)
8 ascii : 头比特数

44 ascii :保留给EDF+
8 ascii : 初始值-1，结束时被赋其他值

8 ascii : 数据持续记录时间，s
4 ascii : 记录几种信号种类

ns * 16 ascii : 电极位置，体温等信息

ns * 80 ascii : 电极信息
ns * 8 ascii : ns * 幅值单位信息

ns * 8 ascii : ns * physical minimum (e.g. -500 or 34)
ns * 8 ascii : ns * physical maximum (e.g. 500 or 40)
ns * 8 ascii : ns * digital minimum (e.g. -2048)
ns * 8 ascii : ns * digital maximum (e.g. 2047)
ns * 80 ascii : 滤波器参数
ns * 8 ascii : 采样率
ns * 32 ascii : 采集信号类型

DATA RECORD
nr of samples[1] * integer : first signal in the data record
nr of samples[2] * integer : second signal
..
nr of samples[ns] * integer : last signal

我们将matlab官方提供的读edf头文件的脚本跑一下：看看是否像如上格式

 [header,~] = edfread('Subject00_1.edf')

“edfRead ” version 2.10 (7.44 KB) by Brett Shoelson：

(温馨提醒：为保证后续阅读，手机用户长按松开拷贝地址到浏览器打开为宜)

http://cn.mathworks.com/matlabcentral/fileexchange/31900-edfread

下载代码，命令行运行结果如下：

我的数据结果：

符合翻译的Header的格式。

帧头之后就是数据。

这里原始数据我用了网上的另一个代码：因为他支持的格式更全，支持edf,rec转mat。

用的是 Alois Schloegl 的脚本。（点击下载源码）

https://files.cnblogs.com/files/myohao/edfsample.zip

mathworks官网也给出了列子：https://ww2.mathworks.cn/matlabcentral/fileexchange/31900-edfread?s_tid=prof_contriblnk

直接运行代码取变量S的前三列，即EDF数据前三通道的脑电数据，画出结果如下：

recorddata = recorddata';
%转置数据,一列为一个通道,采集的数据太长，只展示2000点
plot(recorddata(1:2000,1:3));

将工作区的变量S右键，数据另存为导出为txt。好了有了我们最熟悉的txt格式，接下来选出想分析的通道（列）进行分析就可以了。

（我的数据为一行代表一个通道）

接下来我们跑个简单的算法：FFT分析睡眠数据，我们提取出不同的节律，也就是不同的频段，进行功率谱估计。

以下是四种节律：

α/阿尔法脑波（ALPHA）在大脑中有时出现，有时消失，它并不总是存在。例如，在深睡情况下没有α波；如果一个人在激动状态下，或恐惧，愤怒时，大脑中也没有α脑波。α脑波在初睡或初醒时出现（即半睡半醒时），此时身体处于放松状态，并有自觉的警觉意识。

δ/德尔塔脑波（DELTA）只在深睡时出现。

θ/西塔脑波（THETA）在浅睡时出现。

β/贝塔脑波（BETA）在清醒时出现，伴有需努力能够达到的注意力集中。

不懂傅里叶变换原理的先来这里：最通俗的讲解传送门：

https://zhuanlan.zhihu.com/p/19759362

Matlab下FFT实现代码实例：

https://cn.mathworks.com/help/matlab/ref/fft.html?searchHighlight=FFT&s_tid=doc_srchtitle

由于添加了噪声，fft之后的数据强度接近0.5和1，但不完全相等

%fft部分实验，生成带噪声的数据，fft变换，再反变换回时域数据
%Y = fft(X,n)如果 X 是矩阵，则 fft(X) 将 X 的各列视为向量，并返回每列的傅里叶变换。
%fft的官方解释https://ww2.mathworks.cn/help/matlab/ref/fft.html?searchHighlight=FFT&s_tid=doc_srchtitle
Fs = 500;            % Sampling frequency                    
T = 1/Fs;             % Sampling period       
L = length(recorddata(:,1));  % Length of signal
t = (0:L-1)*T;        % Time vector
%构造一个信号，其中包含幅值为 0.7 的 50 Hz 正弦量和幅值为 1 的 120 Hz 正弦量。
S = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);
%用均值为零、方差为 4 的白噪声扰乱该信号。
X = S + 2*randn(size(t));
%在时域中绘制含噪信号。通过查看信号 X(t) 很难确定频率分量。
figure; plot(1000*t(1:50),X(1:50))
title('Signal Corrupted with Zero-Mean Random Noise');
xlabel('t (milliseconds)');
ylabel('X(t)');
%%
%计算信号的傅里叶变换。
Y = fft(X);
%计算双侧频谱 P2。然后基于 P2 和偶数信号长度 L 计算单侧频谱 P1。
P2 = abs(Y/L);
P1 = P2(1:L/2+1);%计算单侧频率
P1(2:end-1) = 2*P1(2:end-1);
%定义频域 f 并绘制单侧幅值频谱 P1。与预期相符，由于增加了噪声，幅值并不精确等于 0.7 和 1。一般情况下，较长的信号会产生更好的频率近似值。
f = Fs*(0:(L/2))/L;
figure; plot(f,P1) 
title('Single-Sided Amplitude Spectrum of X(t)')
xlabel('f (Hz)')
ylabel('|P1(f)|')

%%
%现在，采用原始的、未破坏信号的傅里叶变换并检索精确幅值 0.7 和 1.0。
Y = fft(S);
P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);

figure; plot(f,P1) 
title('Single-Sided Amplitude Spectrum of S(t)')
xlabel('f (Hz)')
ylabel('|P1(f)|')

有了FFT的函数，在主函数中调用，分割不同频率，分割不同时间，就得到需要的图片了。比如，我想得到10s-150s的θ/西塔脑波（THETA）值。部分源码在这：

http://www.cnblogs.com/myohao/p/8539001.html

测试结果如下