基于小波分析和机器学习的时间序列分析与识别

研究对象：ECG等时间序列信号

方法：小波变换，简单神经网络

首先导入相关模块，需要安装尺度谱模块：pip install scaleogram

和mat4py模块：pip install mat4py

import numpy as npimport pandas as pdimport pywtimport seaborn as snsimport scaleogram as scg import matplotlib.pyplot as pltimport matplotlib.gridspec as GridSpecfrom mat4py import loadmatfrom scipy.fftpack import fft

导入数据，并提取标签

data = loadmat('ECGData.mat')ecg_total = len(data['ECGData']['Data'])ecg_data = []ecg_labels = []for i in range(0, ecg_total):    ecg_data.append(data['ECGData']['Data'][i])    ecg_labels.append(data['ECGData']['Labels'][i])flat_list_ecg_labels = [item for sublist in ecg_labels for item in sublist]

数据可视化

fig = plt.figure(figsize=(12, 6))grid = plt.GridSpec(3, 1, hspace=0.6)arr_signal = fig.add_subplot(grid[0, 0])chg_signal = fig.add_subplot(grid[1, 0])nsr_signal = fig.add_subplot(grid[2, 0])arr_signal.plot(range(0, len(data['ECGData']['Data'][33]), 1), ecg_data[33], color = 'blue')arr_signal.set_xlim(0, 1000)arr_signal.set_title('ARR Signal')chg_signal.plot(range(0, len(data['ECGData']['Data'][100]), 1), ecg_data[100], color = 'red')chg_signal.set_xlim(0, 1000)chg_signal.set_title('CHG Signal')nsr_signal.plot(range(0, len(data['ECGData']['Data'][150]), 1), ecg_data[150], color = 'green')nsr_signal.set_xlim(0, 1000)nsr_signal.set_title('NSR Signal')

进行傅里叶变换

nn = 160signal_length = 1000full_signal_fft_values = np.abs(fft(ecg_data[nn][:signal_length]))x_values_fft = range(0, len(data['ECGData']['Data'][nn]), 1)[:signal_length]fig = plt.figure(figsize=(12, 6))grid = plt.GridSpec(2, 1,hspace=0.6)full_signal = fig.add_subplot(grid[0, 0])fft_comp = fig.add_subplot(grid[1, 0])full_signal.plot(x_values_fft, ecg_data[nn][:signal_length], color = 'green')full_signal.set_xlim(0, 512)full_signal.set_title('Full Signal')fft_comp.plot(x_values_fft, list(full_signal_fft_values), color = 'purple')fft_comp.set_xlim(0, 512)fft_comp.set_ylim(0, 100)fft_comp.set_title('FFT of full signal')

看下所用的Morlet小波啥个样子

axes = scg.plot_wav('morl', figsize=(12,4))

顺便看下小波族大致包含多少小波类

pywt.families(False)

['Haar',

'Daubechies',

'Symlets',

'Coiflets',

'Biorthogonal',

'Reverse biorthogonal',

'Discrete Meyer (FIR Approximation)',

'Gaussian',

'Mexican hat wavelet',

'Morlet wavelet',

'Complex Gaussian wavelets',

'Shannon wavelets',

'Frequency B-Spline wavelets',

'Complex Morlet wavelets']

顺便再给几个小波的波形及相应的频谱

下面进行小波尺度谱变换

# 选择一个默认的小波scg.set_default_wavelet('morl')nn = 33signal_length = 128# 小波变换的尺度范围scales = scg.periods2scales( np.arange(1, signal_length+1) )x_values_wvt_arr = range(0,len(ecg_data[nn]),1)# 绘制信号fig1, ax1 = plt.subplots(1, 1, figsize=(9, 3.5)); ax1.plot(x_values_wvt_arr, ecg_data[nn], linewidth=3, color='blue')ax1.set_xlim(0, signal_length)ax1.set_title('ECG ARR')# 计算小波时间-尺度谱scg.cws(ecg_data[nn][:signal_length], scales=scales, figsize=(10, 4.0), coi = False, ylabel='Period', xlabel='Time', title='ECG_ARR: scaleogram with linear period'); print('Default wavelet function used to compute the transform:', scg.get_default_wavelet(), '(', pywt.ContinuousWavelet(scg.get_default_wavelet()).family_name, ')')

下面开始进行机器学习识别，首先准备机器学习所需要的数据

arr_list = ecg_data[0:95]chf_list = ecg_data[96:125]nsr_list = ecg_data[126:162]arr_split_256 = [np.array_split(arr_list[ii], 256) for ii in range(95)]arr_flatten = [item for sublist in arr_split_256 for item in sublist]chf_split_256 = [np.array_split(chf_list[ii], 256) for ii in range(29)]chf_flatten = [item for sublist in chf_split_256 for item in sublist]nsr_split_256 = [np.array_split(nsr_list[ii], 256) for ii in range(36)]nsr_flatten = [item for sublist in nsr_split_256 for item in sublist]reduce_size = 500full_1500 = (arr_flatten[0:reduce_size] + chf_flatten[0:reduce_size] + nsr_flatten[0:reduce_size])

创建数据集

#机器学习相干模块from sklearn import preprocessingfrom sklearn.model_selection import train_test_splitfs = len(full_1500[0]) sgn_length = 2000 #信号长度size_dataset = len(full_1500)scales = range(1, fs)waveletname = 'morl' #小波 X_full = np.ndarray(shape=(size_dataset, fs-1, fs-1, 3)) #开始生成数据 for i in range(0, size_dataset): if i % 500 == 0: print (i, 'done!') for j in range(0, 3): signal = full_1500[i] coeff, freq = pywt.cwt(signal, scales, waveletname, 1) X_full[i, :, :, j] = coeff[:,:fs-1]#相应的标签list_ecg_labels_arr = ['ARR']*reduce_size list_ecg_labels_chf = ['CHF']*reduce_size list_ecg_labels_nsr = ['NSR']*reduce_size list_ecg_labels = (list_ecg_labels_arr + list_ecg_labels_chf + list_ecg_labels_nsr)le = preprocessing.LabelEncoder()ecg_labels_encoded = le.fit_transform(list_ecg_labels)X_train, X_test, y_train, y_test = train_test_split(X_full, ecg_labels_encoded, test_size=0.25, random_state=42)

使用尺度谱训练ECG分类器

import sysfrom tensorflow import keras

查看一下训练数据

定义基本的神经网络

num_filter, num_classes = 3, 3model = keras.models.Sequential([ keras.layers.Flatten(input_shape=[fs-1, fs-1, num_filter]), keras.layers.Dense(300, activation='relu'), keras.layers.Dense(100, activation='relu'), keras.layers.Dense(num_classes, activation='softmax')])model.compile(loss='sparse_categorical_crossentropy', optimizer='sgd', metrics=['accuracy'])model.summary()

训练模型

history = model.fit(X_train, y_train, epochs=10, validation_data=(X_test, y_test))model.evaluate(X_test, y_test)def plot_learning_curves(history):    pd.DataFrame(history.history).plot(figsize=(8, 5))    plt.grid(True)    plt.gca().set_ylim(0, 1)    plt.xlabel('epochs', size = 16)    plt.xticks(fontsize=14)    plt.yticks(fontsize=14)    plt.savefig('Figure_epochs.png')       plt.close()plot_learning_curves(history)prediction = model.predict(X_test)pred_classes = prediction.argmax(axis=1)

绘制混淆矩阵

cm = confusion_matrix(y_test, pred_classes)cm_norm = cm.astype('float') / cm.sum(axis=1)[:, np.newaxis]df_cm = pd.DataFrame(cm_norm, ['ARR', 'CHF', 'NSR'], ['ARR', 'CHF', 'NSR'])plt.figure(figsize = (10,6))conf = sns.heatmap(df_cm, annot=True, square=True, annot_kws={'size': 12})conf.set_xlabel('Prediction')conf.set_ylabel('True')

基于小波分析和机器学习的时间序列分析与识别 - 哥廷根数学学派的文章 - 知乎
https://zhuanlan.zhihu.com/p/554097033