肌電信號 聚類 Matlab

• 前言
• 一、步驟
• • 1.數(shù)據(jù)采集
  • 2.特征提取
  • 3. 數(shù)據(jù)降維
  • 3. K-means 聚類
  • 4. 評價指標(DB值）
  • 5. 效果
• 二、數(shù)據(jù)說明

Author:溪云楓
2021.10.28

---

前言

本文簡述生物信號（肌電，腦電）數(shù)據(jù)分析步驟，使用PCA, K-means等方法實現(xiàn)信號聚類。

---

一、步驟

1.數(shù)據(jù)采集

肌電，腦電信號可采用各種生物傳感器進行采集，市面上有很多種，有的還可以直接生成濾波后的數(shù)字信號。例如神念公司的腦電采集模塊等。

2.特征提取

肌電信號特征多樣，可根據(jù)分類效果進行擇優(yōu)提取，本文提取了絕對均值、方差、過零點數(shù)等五個特征。常用特征見下圖：

目前較好的特征有表中的（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的行數(shù)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;endmav_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;endend for r=1:size(L,2)data_test=L(:,r);Z=0;for j=1:length(data_test)-1if data_test(j)*data_test(j+1)<0Z=Z+1;endendZC(r)=Z; end %ZC=ZC*2/length(data_test); %斜率變化數(shù)（SSC） [R, C] = size(L); answ = zeros(1, C); slopes = diff(L); threshold=0.03; for i =1:Ccount = 0;for j = 1: R-2if((((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;endendansw(1,i) = count;SSC=answ; end %SSC=SSC*3/length(data_test); data=[mav_data_;rms_data_;var_data_;ZC;SSC]; % 原始數(shù)據(jù)集 data_kinds=zeros(4,30); for j=1:5n=1;for i=1:256:size(data,2)data_kinds(j,n)=mean(data(j,i:i+255));n=n+1;end end

3. 數(shù)據(jù)降維

將高維數(shù)據(jù)投影到低維，增加分類可靠性。本文使用PCA降維方法。另外，t-SNE方法能提高降維后數(shù)據(jù)可視效果，也常使用。

%特征矩陣歸一化 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'); %求樣本的協(xié)方差矩陣（散步矩陣除以(n-1)即為協(xié)方差矩陣） [V,D]=eigs(covX); %求協(xié)方差矩陣的特征值D和特征向量V meanX=mean(data_kinds'); %樣本均值m %所有樣本X減去樣本均值m，再乘以協(xié)方差矩陣（散步矩陣）的特征向量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; % 模式樣本維數(shù) disp('設(shè)置 5 個聚類中心') k = 5; % 設(shè)有k個聚類中心 PM= pcaData2;% 模式樣本矩陣 N = size(PM,1); figure(); subplot(1,2,1); for(i=1:N)plot(PM(i,1),PM(i,2), '*r'); % 繪出原始的數(shù)據(jù)點hold on end xlabel('X'); ylabel('Y'); title('聚類之前的數(shù)據(jù)點');CC = zeros(k,dim); % 聚類中心矩陣，CC(i,:)初始值為i號樣本向量 D = zeros(N,k); % D(i,j)是樣本i和聚類中心j的距離C = cell(1,k); %% 聚類矩陣，對應(yīng)聚類包含的樣本。初始狀況下，聚類i(i<k)的樣本集合為[i],聚類k的樣本集合為[k,k+1,...N] for i = 1:k-1C{i} = [i]; end C{k} = k:N;B = 1:N; % 上次迭代中，樣本屬于哪一聚類，設(shè)初值為1 B(k:N) = k;for i = 1:kCC(i,:) = PM(i,:); endwhile 1 change = 0;%用來標記分類結(jié)果是否變化% 對每一個樣本i，計算到k個聚類中心的距離for i = 1:Nfor 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);endt = 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 endif change == 0 %分類結(jié)果無變化，則迭代停止break;end% 重新計算聚類中心矩陣CCfor i = 1:kCC(i,:) = 0;iclu = C{i};for j = 1:length(iclu)CC(i,:) = PM( iclu(j),: )+CC(i,:);endCC(i,:) = CC(i,:)/length(iclu);end endsubplot(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 onelseif(B(1,i)==2)plot(PM(i,1),PM(i,2), '*r'); %作出第二類點的圖形hold onelseif(B(1,i)==3)plot(PM(i,1),PM(i,2),'*g'); %作出第三類點的圖形hold onelseif(B(1,i)==4) plot(PM(i,1),PM(i,2),'*y'); %作出第四類點的圖形hold onelseplot(PM(i,1),PM(i,2), '*m'); %作出第五類點的圖形hold onend end xlabel('X'); ylabel('Y'); title('聚類之后的數(shù)據(jù)點');% 打印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); %%%保存每次聚類結(jié)果 end eva = evalclusters(pcaData2,cluster,'DaviesBouldin'); figure(2) plot(eva) title('不同簇數(shù)的評價') legend('指標越小表明聚類效果越好') disp('簇數(shù)為4或5時聚類效果較好。簇數(shù)為6時雖然DB值很小，但有的簇僅有一個樣本，聚類不準確。')

5. 效果

說明：簇數(shù)為6時，分類結(jié)果顯示有的類只有一個數(shù)據(jù)點，由于K-means是無監(jiān)督學習，故簇數(shù)為6不可靠。

二、數(shù)據(jù)說明

肌電采集設(shè)備采樣率為2048Hz，貼在受試者手臂上的通道為256個。數(shù)據(jù)文件中，包含了受試者的N個手勢動作,每個動作截取了1秒的256通道采集的肌電信號，因此對應(yīng)于一個手勢動作，數(shù)據(jù)的形式為2048*256的矩陣。讀者可根據(jù)自己的數(shù)據(jù)自行設(shè)計代碼。

總結(jié)

以上是生活随笔為你收集整理的肌电信号 聚类 Matlab的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。