医学信号处理第二次实验报告

电子科技大学生命科学与技术学院标准实验报告

（实验）课程名称生物医学信号处理

2018-2019-第2学期

电子科技大学教务处制表

一、实验室名称：品学楼B302

二、实验名称：随机信号多角度认知和脑电信号特征的认知

三、实验学时：

四、实验原理：

(一)

1、均值公式为:

方差公式为:

2、一个平稳随机信号中两个时间点上的自相关和自协方差公式为：

3、随机信号的功率谱密度函数定义为自相关函数的傅里叶变换：

（二）

1、傅里叶变换对：频域上与窗函数卷积使得频率谱更加平滑

2、脑电波是一些自发的有节律的神经电活动，其频率变动范围在每秒1－30次之间的，可划分为四个波段：

δ（1－3Hz）：为人的深度睡眠阶段的脑电波。

θ（4－7Hz）：为人的睡眠的初期阶段。

α（8－13Hz）：大脑处于完全放松的精神状态下，或是在心神专注的时候出现的脑电波。

β（14－30Hz）这种脑电波反映的是人类在一种通常的、日常的清醒状态下的脑电波情况。

五、实验目的：

1、了解randn和rand产生序列的特征和区别，能够计算各种样本

数字特征及理解一阶统计量与二阶统计量间的区别。

2、利用周期图法估计信号功率谱，知道相关函数结果与功率谱的

DFT关系。

3、探讨不同窗函数对周期图法信号处理结果的影响。通过数据差

别分析信号区别。

六、实验内容：

（一）随机信号的时域和频域认识

1、在波形产生函数中选randn和rand两种波形发生器，各产生一

段随机信号，请观察它们是什么样的信号，描述它们的时域特征。

2、编制一个程序，计算这两个随机信号的样本数字特征，包括均

值、方差、相关函数（xcorr）、协方差函数（xcov），比较并描述这两个信号一阶和二阶统计量的区别。

3、对以上信号样本计算频数直方图（hist）并估计这两个随机信

号的概率密度函数(ksdensity)及估计他们的概率分布函数

(ksdensity) 。

4、利用周期图法估计这两个信号功率谱，比较并描述它们的频域

特征。

5、查看相关函数结果和功率谱之间是否是一对DFT。

（二）开闭眼脑电信号特征的认识

1、按照（学号后两位*班号）选择一路脑电信号，观察和描述开

眼和闭眼脑电信号的时域波形特征。（数据文件“eegclose.mat”

and “eegopen.mat”，Fs=250 Hz，幅度单位：微伏）

2、使用周期图法对开眼和闭眼的脑电信号进行分析，探讨不同窗

函数对分析结果的影响（矩形窗之外选三种窗）；

3、将某一种窗函数下的开眼和闭眼功率谱图进行比较，找出开眼

与闭眼功率谱上存在的差异。（例如测量脑电delta，theta、alpha，beta四个波段内的功率峰值，采用表格方式列出，这样可以比较开眼和闭眼的功率谱分布的差异）。

4、给出一段文字总结开眼和闭眼脑电信号之间的差异总结。

七、实验器材（设备、元器件）：matlab2014b

八、实验步骤:首先整理思路，写程序并记录运行结果，分析数据并写报告。

九、实验数据及结果分析：

（一）程序：

clear all;

clc;

N=10000;

s_randn = randn(1,N);

s_rand = rand(1,N);

%均值方差

mean_randn = mean(s_randn) mean_rand = mean(s_rand)

V_randn = var(s_randn)

V_rand = var(s_rand)

%相关函数

xg_randn = xcorr(s_randn,'biased'); xg_rand = xcorr(s_rand,'biased');

figure(1)

subplot(3,2,1)

plot(1:N,s_randn)

title('正态分布随机数')

subplot(3,2,2)

plot(1:N,s_rand)

title('均匀分布随机数')

subplot(3,2,3)

plot(1:length(xg_randn),xg_randn) title('正态分布随机数自相关')

subplot(3,2,4)

plot(1:length(xg_rand),xg_rand) title('均匀分布随机数自相关')

xfc_randn = xcov(s_randn,'biased'); xfc_rand = xcov(s_rand,'biased'); subplot(3,2,5)

plot(1:length(xfc_randn),xfc_randn) title('正态分布随机数协方差') subplot(3,2,6)

plot(1:length(xfc_rand),xfc_rand) title('均匀分布随机数协方差')

figure(2)

subplot(3,2,1)

hist (s_randn,100);

title('正态分布随机数频率直方图') subplot(3,2,2)

hist (s_rand,100);

title('均匀分布随机数频率直方图') [f1,s1] = ksdensity(s_randn);

[f2,s2] = ksdensity(s_rand); subplot(3,2,3)

plot(s1,f1)

title('正态分布随机数概率密度曲线') subplot(3,2,4)

plot(s2,f2)

title('均匀分布随机数概率密度曲线')

F1 = zeros(1,length(f1));

F2 = zeros(1,length(f2));

for i = 1:length(f1)

F1(i) = sum(f1(1:i));

end

for i = 1:length(f2)

F2(i) = sum(f2(1:i));

end

subplot(3,2,5)

plot(1:length(F1),F1)

title('正态分布随机数概率分布函数曲线') subplot(3,2,6)

plot(1:length(F2),F2)

title('正态分布随机数概率分布函数曲线') pow_n=abs(fft(s_randn,2*N-1)).^2/N;

xgn_fft=abs(fft(flip(fftshift(xg_randn)))); d1=pow_n-xgn_fft;

pow=abs(fft(s_rand,2*N-1)).^2/N;

xg_fft=abs(fft(flip(fftshift(xg_rand))));

d2=pow-xg_fft;

figure(3)

suptitle('验证相关与功率谱是一对DFT')

subplot(1,2,1)

plot(1:length(pow_n),pow_n,'*',1:length(xgn_fft),xgn_fft,'o') legend('正态序列功率谱','正态序列相关')

subplot(1,2,2)

plot(1:length(pow),pow,'*',1:length(xg_fft),xg_fft,'o') legend('平均序列功率谱','正态序列相关')

disp('正态随机信号的计算误差最大值为:')

max(d1)

disp('均匀随机信号的计算误差最大值为:')

max(d2)

结果：

1、2 ：两个随机信号的样本数字特征

表1 信号的均值及方差表

图1 信号的时域图像和二阶统计量图像

3：频数直方图（hist）并估计这两个随机信号的概率密度函数(ksdensity)及估计他们的概率分布函数(ksdensity) 。

图2 频数直方图概率密度函数及概率分布函数4：相关函数结果和功率谱比较

图3 信号的功率谱和相关函数

正态随机信号的计算误差最大值为: 4.4409e-15

均匀随机信号的计算误差最大值为: 1.3088e-13

（一）程序：

clear all

close=load('eegclose');

open=load('eegopen');

close=close.eegclose(:,8*2);

open=open.eegopen(:,8*2);

n = length(close);

fs = 250;

x=(0:n-1)'/fs;

figure(1)

subplot(2,2,1)

plot(x,close);

xlabel('时间/s');ylabel('幅度/mV');title('闭眼信号') subplot(2,2,2)

plot(x,open);

xlabel('时间/s');ylabel('幅度/mV');title('睁眼眼信号')

%构造

win=cell(1,4);

rect=boxcar(n);win{1,1}=rect;

trian=triang(n);win{1,2}=trian;

hamm=hamming(n);win{1,3}=hamm;

bman=blackman(n);;win{1,4}=bman;

figure(2)

plot(1:n,rect);hold on

plot(1:n,trian);hold on

plot(1:n,hamm);hold on

plot(1:n,bman);

legend('矩形窗','三角窗','海明窗','布莱克曼窗')

close_w = cell(1,4);

open_w = cell(1,4);

for i = 1:length(win)

close_w{1,i }= abs(fft(close.*win{1,i})).^2/n;

[pks1,locs1] = findpeaks(close_w{1,i});

open_w{1,i} = abs(fft(open.*win{1,i})).^2/n;

[pks2,locs2] = findpeaks(open_w{1,i});

end

freq = (0:n/2-1)/n *fs;

figure(3)

names={'矩形窗','三角窗','海明窗','凯撒窗'};

%在功率谱中标注出δ(1-3Hz)、θ(4-7Hz)、α(8-13Hz)、β(14-30Hz) locate_mid = [1.5,5.5,10.5,22];

%生成四种脑电波的定位矩阵

loc = {[1,3] [4,7] [8,13] [14,30]};

for i = 1:4

subplot(4,2,2*i-1)

plot(freq,close_w{1,i}(1:n/2))

xlabel('频率/Hz')

title(['闭眼加',names{1,i},'的功率谱']);

axis([0 40 0 max(close_w{1,i})/10])

set(gca,'XTickmode','manual','Xtick',locate_mid) subplot(4,2,2*i)

plot(freq,open_w{1,i}(1:n/2))

xlabel('频率/Hz')

title(['睁眼加',names{1,i},'的功率谱']);

axis([0 40 0 max(open_w{1,i})/10])

set(gca,'XTickmode','manual','Xtick',locate_mid) end

peaks_c = zeros(4,4);

peaks_o = zeros(4,4);

locate = cell(1,4);

for j=1:4

for i = 1:4

locate{1,i} = find ((freq > loc{1,i}(1)) & (freq < loc{1,i}(2)));

peaks_c(j,i) = max( findpeaks( close_w{1,j}( (locate{1,i}) ) ) );

peaks_o(j,i) = max( findpeaks( open_w{1,j}( (locate{1,i}) ) ) );

end

end

f = figure('Position',[440 500 461 146]);

% Create the column and row names in cell arrays

cnames = {'δ(1-3Hz)','θ(4-7Hz)','α(8-13Hz)','β(14-30Hz)'};

% Create the uitable

t = uitable(f,'Data',peaks_c, 'ColumnName',cnames,'RowName',names);

% Set width and height

t.Position(3) = t.Extent(3);

t.Position(4) = t.Extent(4);

f = figure('Position',[440 500 461 146]);

% Create the column and row names in cell arrays

cnames = {'δ(1-3Hz)','θ(4-7Hz)','α(8-13Hz)','β(14-30Hz)'};

% Create the uitable

t = uitable(f,'Data',peaks_o, 'ColumnName',cnames,'RowName',names); t.Position(3) = t.Extent(3);

t.Position(4) = t.Extent(4);

结果：

1：开眼和闭眼脑电信号的时域波形特征

图4 闭眼和开眼信号时域图像

2：使用周期图法对开眼和闭眼的脑电信号进行分析

图5闭眼和睁眼信号加窗分析信号功率谱

3、比较开眼和闭眼的功率谱分布的差异

表2 闭眼信号各各波段功率谱峰值

表3 睁眼信号各各波段功率谱峰值

表4 两信号各各波段功率谱峰值差值表

图7闭眼信号各波段峰值对比图

矩形窗

三角窗 海明窗 布莱克曼窗 0.00E+00

1.00E+03

2.00E+03

3.00E+03

4.00E+03

5.00E+03 δ(1-3Hz)

θ(4-7Hz)

α(8-13Hz)

β(14-30Hz)

闭眼信号各波段峰值

矩形窗 三角窗

海明窗

布莱克曼窗

图8 闭眼信号各波段峰值对比图

图9 闭、睁眼信号各波段峰值差值对比图

（闭眼减去睁眼）

十、总结及心得体会： （一）

1、正态分布随机信号：值为-4到4之间的随机信号，但大多数值

矩形窗

三角窗 海明窗 布莱克曼窗 0.00E+00

1.00E+03

2.00E+03

3.00E+03

4.00E+03

5.00E+03

6.00E+03 δ(1-3Hz)

θ(4-7Hz)

α(8-13Hz)

β(14-30Hz)

睁眼信号各波段峰值

矩形窗 三角窗

海明窗

布莱克曼窗

集中在-2到2之间。可以看到有均匀分布的较为凸出的“毛刺”。

平均分布随机信号：值均为正数，在0到1之间，没有明显凸出的“毛刺”

2、一维统计量：

正态信号均值接近0，方差接近1。

均匀信号均值接近0.5，方差接近0。

二维统计量：

正态信号自相关与自协方差类似，除了中点为1外值越靠近中间越大，但接近于零。

均匀信号自相关仍是中点凸出，其他部分形状从两边向中间线性递增形成类似于三角形形状。自协方差形状与正态分布类似。3、正态信号：频数分布成正态分布特征，概率密度函数(ksdensity)形状类似。概率分布函数变化集中在中位数附近。

均匀信号：频数分布成均匀分布特征，概率密度函数(ksdensity)形状类似，但出现了小于0和大于1的部分。概率分布函数变化均匀几乎呈线性。

4、正态信号：功率谱上的点频率分布（横向分布）平均，纵向分布功率越小的点越多。

均匀信号：功率功率谱基本上全为0。

5、计算的差值的峰值在10E-11以下，从图像也可以观察到相关函数结果和功率谱之间是一对DFT。（由于相关函数不完全对称，很小的虚部忽略）

（二）

1、闭眼信号：时间为18s，幅度在-40mv到30mv之间。毛刺幅

度较睁眼信号密，且幅度略大。

睁眼信号：时间18s，幅度在-30mv到40mv之间。可以观测到第6s的时候又一次剧烈变化。

2、通过图7和图8各波段峰值比较图可以明显看出，窗函数会使

得频率谱变得平滑。可以减少泄露但也会降低频谱分辨率。

3、通过波段峰值差值比较图图9和差值表4可以明显看出：

睁眼和闭眼波形的差别主要在δ(1-3Hz) ，α(8-13Hz)。如果将差值进行归一化，取差值与峰值平均的比值为纵坐标，那么差别最大的将会是α波。

这也许说明被试闭眼时处于静息状态，甚至快要睡着了。

十一、对本实验过程及方法、手段的改进建议：无