Python 替代 Matlab 做数据分析——数据重采样与数据分割

数据的预处理在数据分析中占有非常重要的比例，本文数据来源和 Matlab 代码来自台 湾阳明交通大学的卢家锋老师的Matlab 图形使用者界应用于生医讯号分析。

本文学习如何分别用 Matlab 和 Python 对时序信号进行基础的操作，包括加载数据，数据重采样(Resampling)，分割数据，平均数据以及绘图等。

1、数据简介

本文中用到了两个数据，其中EMG_ICA.mat 为EMG(肌电信号)，另一个fNIRSdata.mat为EEG(脑电)信号。第一个项目是读取EMG信号，并对其进行重采样。

在第二个项目中，如下图所示，受试者根据指令在特定时间点抬起左臂或者右臂，并记录EEG信号。该项目的目的是提取抬手臂的时间戳，分割信号。

2、数据加载

对于 .mat 的数据文件，Matlab 直接使用load()函数直接加载就可以了，加载之后相应的变量就自动加载到内存中了。在 python 中，需要使用 scipy.io.loadmat()函数，数据会以dict 的格式加载

# Matlab
load('EMG_ICA.mat')  # the EMG data are saved as "fdata" matrix

# Python
mat_data = scipy.io.loadmat('EMG_ICA.mat')
data = mat_data['fdata'].T

3、 数据重采样

在 Matlab 中使用 resample(org_signal,p,q)对数据的进行重采样，其中p,q 为新旧采样率的整数比例，可以使用 [p,q]=rat(new_SR1/org_SR)函数获得。在 Python 中，需要使用scipy.signal.resample(x, num) 其中 x 为原始信号，num 为重采样的采样数量，可以使用下面这个公式计算获得

num_new_samples=num_org_samples*new_sampling_rate/org_sample_rate

时间轴可以由以下公式计算：

# Matlab
org_taxis=[1:length(fdata)]/org_SR
# Python
org_taxis = np.arange(len(org_signal))/org_SR

# Matlab
[p,q]=rat(new_SR1/org_SR); # get the two integer matrices that p/q aproximate new_SR1/org_SR
new_signal1=resample(org_signal,p,q);

# Python
org_signal = data[0]
num_samples1 = len(org_signal)*new_SR1//org_SR
new_signal1 = signal.resample(org_signal, num_samples1)

4、找切割点

在 Matlab 中可以使用 find() 配合括号内的逻辑运算来寻找切割点。在 Python 中需要使用np.nonzero()函数。

注意：
find() 和 np.nonzero() 只能获取 index。
np.nonzero() 返回的是一个 tuple, 需要索引[0] 获取 index

# Matlab
right_time=find(stimarker(:,1)==1);
left_time=find(stimarker(:,2)==1);

# Phython
right_time = np.nonzero(stimarker[:,0] == 1)[0]
left_time = np.nonzero(stimarker[:,1] == 1)[0]

5、 画图

Matlab 中直接用 plot 函数即可，如果需要在同一张 figure 上继续画图，需要 ·hold on·，Python 中只要不申明新的 plt.figure() 就可以在同一张 figure 中继续画图。

如果需要在 (0,x) 上画竖线从y_min到y_max, 仅需要
plot(x = [x, x], y=[y_min,y_max])
同理，在(0,y)上画横线
plot(x=[x_min, x_max], y=[y,y])

# Matlab
figure, plot(taxis,signal), hold on
for i=1:length(right_time)
   plot([taxis(right_time(i)) taxis(right_time(i))],[min(signal) max(signal)],'r--','linewidth',2) % red for right
end
for i=1:length(left_time)
   plot([taxis(left_time(i)) taxis(left_time(i))],[min(signal) max(signal)],'k--','linewidth',2) % black for left
end
xlabel('time(s)')
ylabel('HbO2 concentration (a.u.)')

# Phython
plt.figure()
plt.plot(taxis, signal)
for i in range(right_time.size):
    plt.plot([taxis[right_time[i]],taxis[right_time[i]]], [np.min(signal),np.max(signal)],'r--')

for i in range(left_time.size):
    plt.plot([taxis[left_time[i]],taxis[left_time[i]]], [np.min(signal),np.max(signal)],'k--')

plt.xlabel('time(s)')
plt.ylabel('HbO2 concentration(a.u.)')

如下图所示，抬左手的时刻和抬右手的时刻分别用红色和黑色在EEG信号中标注了出来。

6、切割并分类

最后一步，我们将信号按照上图的标记进行切割，并台左手和右手之后对应的 ECG 信号分开展示。

# Matlab
right_Hbsegment=zeros(segment_tplength,length(right_time));
for i=1:length(right_time)
    right_Hbsegment(:,i)=signal(right_time(i):right_time(i)+segment_tplength-1);
end

left_Hbsegment=zeros(segment_tplength,length(left_time));
for i=1:length(left_time)
    left_Hbsegment(:,i)=signal(left_time(i):left_time(i)+segment_tplength-1);
end

# Python
right_Hbsegment=np.zeros([segment_tplength, right_time.size])
for i in range(right_time.size):
    right_Hbsegment[:,i] = signal[right_time[i]:right_time[i]+segment_tplength]

left_Hbsegment=np.zeros([segment_tplength, left_time.size])
for i in range(left_time.size):
    left_Hbsegment[:,i] = signal[left_time[i]:left_time[i]+segment_tplength]

如下图所示，抬左手和右手所引起的脑电变化如下。

为了更好地分析，抬左右手对 EEG 的不同影响，我们将5组数据求取平均数并放在同一张图中进行比较。

# Matlab
right_BlockAvg=mean(right_Hbsegment,2);
left_BlockAvg=mean(left_Hbsegment,2);

figure,
plot(taxis(1:segment_tplength),right_BlockAvg,'b'),hold on
plot(taxis(1:segment_tplength),left_BlockAvg,'r'),hold on

xlabel('time(s)')
ylabel('HbO2 concentration (a.u.)')
title('Arm lifting')
legend('right-arm','left-arm')

# Python
right_BlockAvg=np.mean(right_Hbsegment,1)
left_BlockAvg=np.mean(left_Hbsegment,1)

plt.figure()
plt.plot(taxis[:segment_tplength],right_BlockAvg, label='right_arm')
plt.plot(taxis[:segment_tplength],left_BlockAvg, label='left_arm')
plt.legend()
plt.xlabel('time(s)')
plt.ylabel('HbO2 concentration (a.u.)')
plt.title('Arm lifting')

如下图所示，经过平均处理之后，抬不同的手对 ECG 的变化是不是有了明显的区别。

文章转载自公众号：Tensorflow机器学习