2021-10-28

2023-10-11 20:44:52

肌电信号聚类 Matlab

前言
一、步骤
二、数据说明

Author:溪云枫
2021.10.28

前言

本文简述生物信号（肌电，脑电）数据分析步骤，使用PCA, K-means等方法实现信号聚类。

一、步骤

1.数据采集

肌电，脑电信号可采用各种生物传感器进行采集，市面上有很多种，有的还可以直接生成滤波后的数字信号。例如神念公司的脑电采集模块等。

2.特征提取

肌电信号特征多样，可根据分类效果进行择优提取，本文提取了绝对均值、方差、过零点数等五个特征。常用特征见下图：

目前较好的特征有表中的（1,4,13,15,18,23,31,33,38,39,48）。

代码如下（Matlab）：

length_t=100; %窗长 100 ms
overlap=50; % overlap 50 ms
for r=1:size(L,2)
    j=1;
    data_test=L(:,r);
    for i=1:overlap:size(data_test,1)   % size(thisdata,1) 返回thisdata的行数
        if i+length_t>size(data_test,1)
            mav_data_(r)=mean(mav_data);
            rms_data_(r)=mean(rms_data);
            var_data_(r)=mean(var_data);
            break;
        end
        mav_data(j)=(sum(abs(data_test(i:i+length_t))))/length_t;
        rms_data(j)=sqrt((sum((data_test(i:i+length_t)).^2))/length_t);
        mean_thisdata=mean(data_test(i:i+length_t));
        var_data(j)=(sum((data_test(i:i+length_t)-mean_thisdata).^2))/length_t;
        j=j+1;
    end
  
end
for r=1:size(L,2)
    data_test=L(:,r);
    Z=0;
    for j=1:length(data_test)-1
        if data_test(j)*data_test(j+1)<0
           Z=Z+1;
        end
    end
    ZC(r)=Z;     
end
%ZC=ZC*2/length(data_test);
%斜率变化数（SSC）
[R, C] = size(L);
answ = zeros(1, C);
slopes = diff(L);
threshold=0.03;
for i =1:C
    count = 0;
    for j = 1: R-2
        if((((slopes(j,i) > 0) && (slopes(j+1,i) < 0))...
                || ((slopes(j,i) < 0) && (slopes(j+1,i) > 0)))...
                && ((abs(slopes(j,i)) >= threshold)...
                || (abs(slopes(j+1,i)) >= threshold)))
        count = count + 1;
        end
    end
    answ(1,i) = count;
    SSC=answ;
end
%SSC=SSC*3/length(data_test);
data=[mav_data_;rms_data_;var_data_;ZC;SSC];  % 原始数据集
data_kinds=zeros(4,30);
for j=1:5
    n=1;
    for i=1:256:size(data,2)
        data_kinds(j,n)=mean(data(j,i:i+255));
        n=n+1;
    end
end

3. 数据降维

将高维数据投影到低维，增加分类可靠性。本文使用PCA降维方法。另外，t-SNE方法能提高降维后数据可视效果，也常使用。

%特征矩阵归一化
data_kinds = data_kinds./repmat(sqrt(sum(data_kinds.^2,1)),size(data_kinds,1),1);

% PCA
f=2;              % 五维降至二维
[Row,Col]=size(data_kinds');
covX=cov(data_kinds');  %求样本的协方差矩阵（散步矩阵除以(n-1)即为协方差矩阵）
[V,D]=eigs(covX);  %求协方差矩阵的特征值D和特征向量V
meanX=mean(data_kinds');         %样本均值m
%所有样本X减去样本均值m，再乘以协方差矩阵（散步矩阵）的特征向量V，即为样本的主成份SCORE
tempX= repmat(meanX,Row,1);
SCORE2=(data_kinds'-tempX)*V;       %主成份：SCORE
pcaData2=SCORE2(:,1:f);   %降维后的二维特征矩阵

3. K-means 聚类

代码：

% 分类
% k-means聚类算法
% main variables
dim = 2; % 模式样本维数
disp('设置 5 个聚类中心')
k = 5;   % 设有k个聚类中心
PM= pcaData2;% 模式样本矩阵 
N = size(PM,1);
figure();
subplot(1,2,1);
for(i=1:N)
   plot(PM(i,1),PM(i,2), '*r'); % 绘出原始的数据点
   hold on
end
xlabel('X');
ylabel('Y');
title('聚类之前的数据点');

CC = zeros(k,dim); % 聚类中心矩阵，CC(i,:)初始值为i号样本向量
D = zeros(N,k); % D(i,j)是样本i和聚类中心j的距离

C = cell(1,k); %% 聚类矩阵，对应聚类包含的样本。初始状况下，聚类i(i<k)的样本集合为[i],聚类k的样本集合为[k,k+1,...N]
for i = 1:k-1
    C{i} = [i];
end
C{k} = k:N;

B = 1:N; % 上次迭代中，样本属于哪一聚类，设初值为1
B(k:N) = k;

for i = 1:k
    CC(i,:) = PM(i,:);
end

while 1   
    change = 0;%用来标记分类结果是否变化
    % 对每一个样本i，计算到k个聚类中心的距离
    for i = 1:N
        for j = 1:k
%             D(i,j) = eulerDis( PM(i,:), CC(j,:) );
              D(i,j) = sqrt((PM(i,1) - CC(j,1))^2 + (PM(i,2) - CC(j,2))^2);
        end
        t = find( D(i,:) == min(D(i,:)) ); % i属于第t类
        if B(i) ~= t % 上次迭代i不属于第t类
            change = 1;
            % 将i从第B(i)类中去掉
            t1 = C{B(i)};
            t2 = find( t1==i );            
            t1(t2) = t1(1);
            t1 = t1(2:length(t1)); 
            C{B(i)} = t1;

            C{t} = [C{t},i]; % 将i加入第t类

            B(i) = t;
        end        
    end

    if change == 0 %分类结果无变化，则迭代停止
        break;
    end

    % 重新计算聚类中心矩阵CC
    for i = 1:k
        CC(i,:) = 0;
        iclu = C{i};
        for j = 1:length(iclu)
            CC(i,:) = PM( iclu(j),: )+CC(i,:);
        end
        CC(i,:) = CC(i,:)/length(iclu);
    end
end

subplot(1,2,2);
plot(CC(:,1),CC(:,2),'o')
hold on
for(i=1:N)
    if(B(1,i)==1)
        plot(PM(i,1),PM(i,2),'*b'); %作出第一类点的图形
        hold on
    elseif(B(1,i)==2)
        plot(PM(i,1),PM(i,2), '*r'); %作出第二类点的图形
     hold on
    elseif(B(1,i)==3)
        plot(PM(i,1),PM(i,2),'*g'); %作出第三类点的图形
        hold on
    elseif(B(1,i)==4)    
        plot(PM(i,1),PM(i,2),'*y'); %作出第四类点的图形
        hold on
    else
        plot(PM(i,1),PM(i,2), '*m'); %作出第五类点的图形
        hold on
    end
end
xlabel('X');
ylabel('Y');
title('聚类之后的数据点');
      % 打印C,CC     
for i = 1:k    %输出每一类的样本点标号
    str=['第' num2str(i) '类包含点:  ' num2str(C{i})];
    disp(str);
end

4. 评价指标(DB值）

cluster = zeros(size(pcaData2,1),6);
for i=1:6
cluster(:,i) = kmeans(pcaData2,i,'replicate',50); %%%保存每次聚类结果
end
eva = evalclusters(pcaData2,cluster,'DaviesBouldin'); 
figure(2)
plot(eva)
title('不同簇数的评价')
legend('指标越小表明聚类效果越好')
disp('簇数为4或5时聚类效果较好。簇数为6时虽然DB值很小，但有的簇仅有一个样本，聚类不准确。')

5. 效果

说明：簇数为6时，分类结果显示有的类只有一个数据点，由于K-means是无监督学习，故簇数为6不可靠。

二、数据说明

肌电采集设备采样率为2048Hz，贴在受试者手臂上的通道为256个。数据文件中，包含了受试者的N个手势动作,每个动作截取了1秒的256通道采集的肌电信号，因此对应于一个手势动作，数据的形式为2048*256的矩阵。读者可根据自己的数据自行设计代码。

码农公寓

肌电信号 聚类 Matlab

前言