python玩轉(zhuǎn)信號(hào)處理與機(jī)器學(xué)習(xí)入門

作者：王鎮(zhèn)

面對毫無規(guī)律的隨機(jī)信號(hào)，看著雜亂無章的振動(dòng)波形，你是否也像曾經(jīng)的我一樣一頭霧水，不知從何處下手。莫慌，接下來小編就帶你入門怎樣用python處理這些看似毫無卵用實(shí)則蘊(yùn)藏巨大信息的隨機(jī)信號(hào)。我們?nèi)粘Ｉ钪兴姷男碾妶D，聲波圖都是信號(hào)在時(shí)域上的一種表現(xiàn)，但它們無法呈現(xiàn)出信號(hào)在頻域上的信息。因此，本文將主要介紹信號(hào)從時(shí)域到頻域上的一些變換，常見的有FFT(快速傅里葉變換)，PSD(功率譜密度)，auto-correlation(自相關(guān)分析)。最后小編將帶你完成一個(gè)實(shí)例，通過手機(jī)采集的振動(dòng)信號(hào)識(shí)別人體的動(dòng)作。

一、介紹

本部分將介紹FFT，PSD，auto-correlation的基本概念以及python代碼實(shí)現(xiàn)。

1.1 混合信號(hào)

圖1 信號(hào)在時(shí)域上的表現(xiàn)

圖2 信號(hào)在頻域上的表現(xiàn)

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

1.2 FFT

FFT英文全稱Fast Fourier Transformation，即快速傅里葉變換，它可以輕松地分析出混合信號(hào)中的各頻率組成成分。對上述中的混合信號(hào)做FFT變換，結(jié)果如圖2(a)，可以明顯地看到混合信號(hào)包含頻率分別為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一樣，反映的是信號(hào)在頻域上的信息。其中PSD頻譜圖脈沖下方的面積表示信號(hào)在該頻率上的能量分布。

對上述中的混合信號(hào)做PSD變換，結(jié)果如圖2(b)，可以明顯地看到混合信號(hào)在頻率為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是自相關(guān)的意思，它可以求出信號(hào)的自相關(guān)性，即信號(hào)經(jīng)過一個(gè)時(shí)延后與自身的相似性。對上述中的混合信號(hào)計(jì)算Autocorrelation，結(jié)果如圖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

二 實(shí)例

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

圖3 數(shù)據(jù)集分布

2.1 數(shù)據(jù)可視化

隨機(jī)選取一信號(hào)，繪出其在時(shí)域和頻域上的波形圖如下所示，繪圖代碼詳見項(xiàng)目鏈接：

圖4 一個(gè)例子的展示

2.2 特征提取

做好了頻譜變換之后，我們需要從中提取特征，這樣才能應(yīng)用我們所熟悉的諸如隨機(jī)森林，邏輯回歸，支持向量機(jī)之類的機(jī)器學(xué)習(xí)模型。那么提取什么特征呢，一種方式是提取脈沖(peak)發(fā)生時(shí)所在的橫縱坐標(biāo)，我們提取頻譜中的前5個(gè)脈沖的橫縱坐標(biāo)作為特征。其中提取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個(gè)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 特征提取詳細(xì)介紹

2.3 模型訓(xùn)練及結(jié)果展示

當(dāng)構(gòu)建完特征矩陣以及其對應(yīng)的標(biāo)簽之后，我們可以選擇scikit-learn庫中的相關(guān)模型進(jìn)行訓(xùn)練。

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))

結(jié)果如下。

圖6 分類結(jié)果展示

圖7 模型比較

正如結(jié)果所展示的那樣，我們能以相當(dāng)高的準(zhǔn)確率(89%)對這些信號(hào)進(jìn)行分類，取得這個(gè)結(jié)果，我們甚至都沒有做任何手動(dòng)的特征工程，所有特征都是自動(dòng)獲取的，對于每個(gè)變換，我們?nèi)∏拔鍌€(gè)峰值的橫縱坐標(biāo)(若少于五個(gè)則填充0)。可以理解的一點(diǎn)是，這270個(gè)特性中的一些特征將比其他特征提供更多的信息，若我們能主動(dòng)地選擇對分類很重要的特征進(jìn)行組合，或者對超參數(shù)進(jìn)行優(yōu)化，相信準(zhǔn)確率還能繼續(xù)提高。

參考文獻(xiàn)：

關(guān)于我們

Mo(網(wǎng)址：https://momodel.cn) 是一個(gè)支持 Python的人工智能在線建模平臺(tái)，能幫助你快速開發(fā)、訓(xùn)練并部署模型。

近期 Mo 也在持續(xù)進(jìn)行機(jī)器學(xué)習(xí)相關(guān)的入門課程和論文分享活動(dòng)，歡迎大家關(guān)注我們的公眾號(hào)獲取最新資訊！

來源：oschina

鏈接：https://my.oschina.net/u/4109778/blog/4277567

總結(jié)

以上是生活随笔為你收集整理的python信号处理教程_python玩转信号处理与机器学习入门的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

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