當我們將肌電訊號轉換成頻譜(Power spectrum)後,我們的確可
以看到頻率和其強度的分佈,不過這代表什麼意義呢?
我們知道我們的肌肉由不同型態組成,簡單的說可以分為快肌
(fast twitch)跟慢肌(slow twitch),快肌肌肉纖維較粗,動作速
度快,不過較容易疲勞,慢肌動作速度慢,但是可以作用持續時間長,
就傳遞肌肉表面膜電位而言,快肌就像比較粗的電線,傳遞訊號的速度
較快,因此由電極片所接收到單位時間內的電位差訊號相對而言就會比
較多,所以頻率會落在比較高的位置,反觀慢肌,則為相反地現象。
而我們的肌肉,則由不同比例的快慢肌所組成,例如我們強而有力
的肱二頭肌(biceps),快肌的組成較多,所以也比較容易酸。
不過因為頻譜看出來的是各個頻率強度的變化,如果要相互比較的
話,則必須要像一個常態分佈圖形一樣,有平均數(mean),等數值,
比較常用來作為計算的為中位頻率(Median frequency)。
中位頻率(Median frequency),其定義為頻譜分析圖中,將其強
度由低頻積分到高頻,積分面積為全部面積一半時的頻率。因此就可以
中位頻率來作為訓練或是疲勞測試後,改變肌肉組成比例(訓練)或是
影響肌肉放電(疲勞)等改變觀察。
因此利用傅立葉轉換以及if指令來求得中位頻率,當積分到某一點
時,大於全部等積分的一半時的頻率即是中位頻率。
由於肌電訊號是時間性的,我們比較常做的方式會像之前所用的以
window的方式,每多少點計算一次,重複多少點,然而取樣頻率
(sampling rate)會影響最後每一個點所代表的頻率,所以也必須考慮
。
function [med_fq]=fft_i(data,n,op,sr)
% fft for n point
% op: the point of overlapping
% med_fq is the median frequency of the power spectrum
% sr is the sampling frequency
% calculate the times of the loop (-1 due to the initial loop is 0)
l=fix(length(data)/op-1);
fq=(1:n/2).*sr/n; % calculate the frequency of each point
for i = 0:(l-1)
int=1+op*i; % the initial value of the loop
f=fft(data(int:(int+511))); % fft for the data
mag=abs(f); % Amplitude of the fft
mag=mag(1:n/2);
for k=1:n/2
if sum(mag(1:k)) >= sum(mag)/2
med=k;
break
end
end
med_fq(i+1,1)=k*sr/n;
end
% calculate the final part of the sample
int=1+op*l;
f=fft(data(int:end)); % fft for the data
mag=abs(f); % Amplitude of the fft
mag=mag(1:n/2);
for k=1:n/2
if sum(mag(1:k)) >= sum(mag)/2
med=k;
break
end
end
med_fq(l+1,1)=k*sr/n;
首先必須先計算迴圈的次數,以所有的點數除上重複的點數(通常
會以window點數的一半)減一(種樹法則),然後就每多少點作一次傅
立葉轉換並且求得中位頻率,值得注意的是,因為常常所選取的資料數
並不能剛好夠最後一筆,舉例而言,也就是說原本應該要以512個點作計
算,但是迴圈到最後一次時可能只剩下500點,由於fft這個函數當不夠
點的時候會比剩下的點以零帶入,所以最後一個窗口則以另外方式計算
。
執行結果
>> data=textread('D:/EMG_fft.txt');
>> plot(data);
>> med_fq=fft_i(data,512,256,1000);
>> figure(2);plot(med_fq)
這是原始資料圖,資料可以由此下載
經由計算後為
一開始很高可能是因為在還沒有任何訊號時就已經存在的高頻雜訊,
此時就要利用之前提得濾波器給予濾波。
0 意見:
張貼留言