作者：王鎮

面對毫無規律的隨機信號，看著雜亂無章的振動波形，你是否也像曾經的我一樣一頭霧水，不知從何處下手。莫慌，接下來小編就帶你入門怎樣用python處理這些看似毫無卵用實則蘊藏巨大信息的隨機信號。我們日常生活中所見的心電圖，聲波圖都是信號在時域上的一種表現，但它們無法呈現出信號在頻域上的信息。因此，本文將主要介紹信號從時域到頻域上的一些變換，常見的有FFT(快速傅里葉變換)，PSD(功率譜密度)，auto-correlation(自相關分析)。最后小編將帶你完成一個實例，通過手機采集的振動信號識別人體的動作。

一、介紹

本部分將介紹FFT，PSD，auto-correlation的基本概念以及python代碼實現。

1.1 混合信號

圖1 信號在時域上的表現

圖2 信號在頻域上的表現

上圖展示了混合信號在時域上的表現形式，圖(a)為一頻率為1Hz，振幅為2的正弦波信號，圖(b)為一頻率為5Hz，振幅為1的正弦波信號，圖(c)為(a)、(b)兩信號的疊加結果。

1.2 FFT

FFT英文全稱Fast Fourier Transformation，即快速傅里葉變換，它可以輕松地分析出混合信號中的各頻率組成成分。對上述中的混合信號做FFT變換，結果如圖2(a)，可以明顯地看到混合信號包含頻率分別為1Hz和5Hz的成分。FFT變換的代碼如下：

from scipy.fftpack import fft

def get_fft_values(y_values, N, f_s):

f_values = np.linspace(0.0, f_s/2.0, N//2)

fft_values_ = fft(y_values)

fft_values = 2.0/N * np.abs(fft_values_[0:N//2])

return f_values, fft_values

1.3 PSD

PSD英文全稱Power Spectral Density，即功率譜密度，它和FFT一樣，反映的是信號在頻域上的信息。其中PSD頻譜圖脈沖下方的面積表示信號在該頻率上的能量分布。

對上述中的混合信號做PSD變換，結果如圖2(b)，可以明顯地看到混合信號在頻率為1Hz和5Hz上的能量分布。PSD變換的代碼如下：

from scipy.signal import welch

def get_psd_values(y_values, N, f_s):

f_values, psd_values = welch(y_values, fs=f_s)

return f_values, psd_values

1.4 Autocorrelation

Autocorrelation是自相關的意思，它可以求出信號的自相關性，即信號經過一個時延后與自身的相似性。對上述中的混合信號計算Autocorrelation，結果如圖2(c)所示

。有趣的是Autocorrelation與PSD是一組FFT變換對，對Autocorrelation作FFT變換可得到PSD，對PSD作IFFT(快速傅里葉逆變換)可得到Autocorrelation。

def autocorr(x):

result = np.correlate(x, x, mode='full')

return result[len(result)//2:]

def get_autocorr_values(y_values, N, f_s):

autocorr_values = autocorr(y_values)

x_values = np.array([ 1.0*jj/f_s for jj in range(0, N)])

return x_values, autocorr_values

二 實例

經過上面的簡單介紹相信你已經了解并掌握了信號在頻域上的變換。寫下來讓我們運用剛學的知識結合機器學習知識來分析一個實例 Human Activity Recognition Using Smartphones Data Set。該數據集通過在30個不同年齡分布的志愿者上做實驗采集得到，在志愿者的腰上固定一手機，手機以50Hz的采樣頻率采集內嵌加速度計和陀螺儀的數據，志愿者被要求做以下六個動作：walking(行走)，walking upstairs(爬樓梯)，walking downstairs(下樓梯)，sitting(坐著)，standing(站著)，laying(躺下)。在濾除噪聲之后，通過滑動窗口的方式將信號切分成2.56s的窗口，每個窗口包含128個樣本。該數據集包含三組數據three-axial linear body acceleration(去除重力加速度的加速度計數據)、three-axial linear total acceleration(包含重力加速度的加速度計數據)、three-axial angular velocity(陀螺儀的數據)，因此共有九個分量，其中訓練集有7392個窗口，測試集有2947個窗口。

圖3 數據集分布

2.1 數據可視化

隨機選取一信號，繪出其在時域和頻域上的波形圖如下所示，繪圖代碼詳見項目鏈接：

圖4 一個例子的展示

2.2 特征提取

做好了頻譜變換之后，我們需要從中提取特征，這樣才能應用我們所熟悉的諸如隨機森林，邏輯回歸，支持向量機之類的機器學習模型。那么提取什么特征呢，一種方式是提取脈沖(peak)發生時所在的橫縱坐標，我們提取頻譜中的前5個脈沖的橫縱坐標作為特征。其中提取peak信息可用detect_peaks。

def get_first_n_peaks(x, y, no_peaks=5):

x_, y_ = list(x), list(y)

if len(x_) >= no_peaks:

return x_[:no_peaks], y_[:no_peaks]

else:#少于5個peaks，以0填充

missing_no_peaks = no_peaks-len(x_)

return x_ + [0]*missing_no_peaks, y_ + [0]*missing_no_peaks

def get_features(x_values, y_values, mph):

indices_peaks = detect_peaks(y_values, mph=mph)

peaks_x, peaks_y = get_first_n_peaks(

x_values[indices_peaks], y_values[indices_peaks])

return peaks_x + peaks_y

def extract_features_labels(dataset, labels, N, f_s, denominator):

percentile = 5

list_of_features = []

list_of_labels = []

for signal_no in range(0, len(dataset)):

features = []

list_of_labels.append(labels[signal_no])

for signal_comp in range(0, dataset.shape[2]):

signal = dataset[signal_no, :, signal_comp]

signal_min = np.nanpercentile(signal, percentile)

signal_max = np.nanpercentile(signal, 100-percentile)

#ijk = (100 - 2*percentile)/10

#set minimum peak height

mph = signal_min + (signal_max - signal_min)/denominator

features += get_features(*get_psd_values(signal, N, f_s), mph)

features += get_features(*get_fft_values(signal, N, f_s), mph)

features += get_features(*get_autocorr_values(signal, N, f_s), mph)

list_of_features.append(features)

return np.array(list_of_features), np.array(list_of_labels)

X_train, Y_train = extract_features_labels(train_signals, train_labels, N, f_s, denominator)

X_test, Y_test = extract_features_labels(test_signals, test_labels, N, f_s, denominator)

圖5 特征提取詳細介紹

2.3 模型訓練及結果展示

當構建完特征矩陣以及其對應的標簽之后，我們可以選擇scikit-learn庫中的相關模型進行訓練。

clf = RandomForestClassifier(n_estimators=1000)

clf.fit(X_train, Y_train)

print("Accuracy on training set is : {}".format(clf.score(X_train, Y_train)))

print("Accuracy on test set is : {}".format(clf.score(X_test, Y_test)))

Y_test_pred = clf.predict(X_test)

print(classification_report(Y_test, Y_test_pred))

結果如下。

圖6 分類結果展示

圖7 模型比較

正如結果所展示的那樣，我們能以相當高的準確率(89%)對這些信號進行分類，取得這個結果，我們甚至都沒有做任何手動的特征工程，所有特征都是自動獲取的，對于每個變換，我們取前五個峰值的橫縱坐標(若少于五個則填充0)。可以理解的一點是，這270個特性中的一些特征將比其他特征提供更多的信息，若我們能主動地選擇對分類很重要的特征進行組合，或者對超參數進行優化，相信準確率還能繼續提高。

參考文獻：

關于我們

Mo(網址：https://momodel.cn) 是一個支持 Python的人工智能在線建模平臺，能幫助你快速開發、訓練并部署模型。

近期 Mo 也在持續進行機器學習相關的入門課程和論文分享活動，歡迎大家關注我們的公眾號獲取最新資訊！

總結

以上是生活随笔為你收集整理的python处理心电图_【技术博客】Python玩转信号处理与机器学习入门的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。