python傅里叶变换 信号处理 序列_基于信号处理技术的机器学习

本文讲解一篇关于基于信号处理技术的机器 学习的博客：《Machine Learning with Signal Processing Techniques》，极力推荐！！

链接：http://ataspinar.com/2018/04/04/machine-learning-with-signal-processing-techniques/ 1、引言

随机信号分析是一个与随机信号的处理、修正和分析相关的学科。

任何有物理或工程背景的人都一定程度上了解信号分析技术，这些技术是什么，以及如何使用它们来分析、建模和分类信号。但是，来自不同领域的数据科学家，比如计算机科学或统计学，可能没有意识到这些技术可以带来的分析能力。

在这篇博客文章中，我们将介绍如何使用随机信号分析技术，结合传统的机器学习分类器，来对时间序列和信号进行精确地建模和分类。

1、信号的基本知识 1.1 信号 VS 时间序列

你可能经常遇到描述数据集的术语，例如时间序列和信号，但是并不清楚它们之间究竟有什么区别。

在时间序列数据集中，待预测值 y 是时间 t 的函数 y=y(t)，但是信号是一个更一般的概念：因变量 y 未必是时间的函数，也可以是空间坐标的函数 y=y(x, y)，或者距原点距离的函数 y=y(r)。信号具有多种多样的形式，像声音信号、图片、雷达数据等。从本质上说，任何东西都可以叫做信号，只要它本身携带信息。

1.2 离散时间信号

连续信号与离散信号的主要区别在于：连续信号有一个连续的自变量(例如，时间或者空间坐标)。如果你将信号放大后观察，你会发现在每一个可以取到实数值的地方都存在连续信号的值。而离散信号只会在某些特定的时刻存在值，比如t=0.1s、t=0.2s、t=0.3s...，在t=0.13s就没有值。

为了能在电脑上分析和可视化一个信号，通常需要将一个模拟连续信号进行数字化(以一定频率进行采样)，将其变为离散信号。

选择一个好的采样率很重要：如果采样率太低，离散信号会丢失很多原始信号的信息。采样率需要满足奈奎斯特采样定律，即采样率要大于原始信号中最大频率的两倍。

2、时域和频域之间的变换(FFT、PSD、Wavelet)

傅里叶分析是一个分析信号周期性的研究领域。如果一个信号含有周期性成分，傅里叶分析可以用来将这个信号分解成周期性成分，并且能够告诉我们每个周期性分量的频率是多少。

傅里叶变换可以将复合信号分解成一系列的周期性成分(三角函数)。

2.1 Python中的快速傅里叶变换

在Python中，一个信号的快速傅里叶变换可以用Scipy库计算。

from scipy.fftpack import fftdef get_fft_values(y_values, T, N, f_s):    f_values = np.linspace(0.0, 1.0/(2.0*T), N//2)    fft_values_ = fft(y_values)    fft_values = 2.0/N * np.abs(fft_values_[0:N//2])    return f_values, fft_valuest_n = 10N = 1000T = t_n / Nf_s = 1/Tx_value = np.linspace(0,t_n,N)amplitudes = [4, 6, 8, 10, 14]frequencies = [6.5, 5, 3, 1.5, 1]y_values = [amplitudes[ii]*np.sin(2*np.pi*frequencies[ii]*x_value) for ii in range(0,len(amplitudes))]composite_y_value = np.sum(y_values, axis=0)f_values, fft_values = get_fft_values(composite_y_value, T, N, f_s)plt.plot(f_values, fft_values, linestyle='-', color='blue')plt.xlabel('Frequency [Hz]', fontsize=16)plt.ylabel('Amplitude', fontsize=16)plt.title("Frequency domain of the signal", fontsize=16)plt.show()

该复合信号是由频率分别为6.5、5、3、1.5和1Hz以及相应振幅为4、6、8、10和14的sin函数组成的。

2.2 Python中的PSD

与傅里叶变换密切相关的一个概念是功率谱密度分析。和傅里叶变换类似，它也是描述一个信号频谱的方法，并且考虑了每个频率的功率分布。一般来说，PSD图中峰值对应的频率和FFT图中峰值对应的频率一样，但是峰值的高度和宽度不一样。PSD图中峰值下面的面积对应着那个频率的功率分布。

from scipy.signal import welchdef get_psd_values(y_values, T, N, f_s):    f_values, psd_values = welch(y_values, fs=f_s)    return f_values, psd_valuest_n = 10N = 1000T = t_n / Nf_s = 1/Tx_value = np.linspace(0,t_n,N)amplitudes = [4, 6, 8, 10, 14]frequencies = [6.5, 5, 3, 1.5, 1]y_values = [amplitudes[ii]*np.sin(2*np.pi*frequencies[ii]*x_value) for ii in range(0,len(amplitudes))]composite_y_value = np.sum(y_values, axis=0)f_values, psd_values = get_psd_values(composite_y_value, T, N, f_s)plt.plot(f_values, psd_values, linestyle='-', color='blue')plt.xlabel('Frequency [Hz]')plt.ylabel('PSD [V**2 / Hz]')plt.show()

2.3 在Python中计算自相关系数

自相关函数计算信号与自身延时版本的相关性：如果一个信号的周期为t，那么该信号和其自身延迟t秒的版本之间就会有很高的相关性。

Scipy库中没有现成的计算自相关系数的函数，但是我们可以借助numpy库中的correlate()函数实现。自相关系数是与时延t有关的函数值，因此t的取值不能超过原始信号的长度。

def autocorr(x):    result = np.correlate(x, x, mode='full')    return result[len(result)//2:]def get_autocorr_values(y_values, T, N, f_s):    autocorr_values = autocorr(y_values)    x_values = np.array([T * jj for jj in range(0, N)])    return x_values, autocorr_valuest_n = 10N = 1000T = t_n / Nf_s = 1/Tx_value = np.linspace(0,t_n,N)amplitudes = [4, 6, 8, 10, 14]frequencies = [6.5, 5, 3, 1.5, 1]y_values = [amplitudes[ii]*np.sin(2*np.pi*frequencies[ii]*x_value) for ii in range(0,len(amplitudes))]composite_y_value = np.sum(y_values, axis=0)t_values, autocorr_values = get_autocorr_values(composite_y_value, T, N, f_s)plt.plot(t_values, autocorr_values, linestyle='-', color='blue')plt.xlabel('time delay [s]')plt.ylabel('Autocorrelation amplitude')plt.show()

通过对自相关系数进行傅里叶变换，可以得到与FFT一样的频率成分。自相关系数与PSD也是一对傅里叶变换，也就是说通过对自相关系数进行傅里叶变换，可以得到PSD；对PSD做傅里叶反变换，可以得到自相关系数。

2.4 小波变换

小波变换适用于动态频谱的信号(信号的频率随着时间不断变换)。

https://www.mathworks.com/help/wavelet/examples.html这个链接里提供了matlab小波工具箱里的一些实用范例，包含了如下内容：

3、统计参数估计和特征提取

前文已经讲述了利用FFT、PSD和自相关函数三种方法来计算信号的频域特性。这三种方法将时域信号转换到频域上，并给出了频谱图。将信号转换到频域上后，我们可以从转换后的信号中提取特征，并将其作为分类器的输入。 那么，可以提取哪些特征呢？一个比较好的想法是将这三种方法得到的图中峰值的坐标作为特征，如下图所示。

3.1 举例：对人类活动数据集进行分类

关于数据集的描述：This dataset contains measurements done by 30 people between the ages of 19 to 48. The measurements are done with a smartphone placed on the waist while doing one of the following six activities:    walking,    walking upstairs,    walking downstairs,    sitting,    standing or    laying.    The measurements are done at a constant rate of 50 Hz. After filtering out the noise, the signals are cut in fixed-width windows of 2.56 sec with an overlap of 1.28 sec.Each signal will therefore have 50 x 2.56 = 128 samples in total.The smartphone measures three-axial linear body acceleration, three-axial linear total acceleration and three-axial angular velocity. So per measurement, the total signal consists of nine components

划分好训练集和测试集后，训练集的大小为(7352，128，9)，测试集的大小为(2497，128，9)。也就是说，训练集中有7352个信号，每个信号的9个分量均含有128个样本点。选择其中的一个信号，对其分别做FFT、PSD和自相关系数变换后的结果如下图所示。

3.2 对训练集和测试集中的所有信号提取特征

特征提取流程：    transform a signal by means of the FFT, PSD or autocorrelation function.    locate the peaks in the transformation with the peak-finding function.所提取的特征个数决定于下面几个因素：· The number of columns in each matrix depend on depends on your choice of features. · Each signal has nine components, and for each component you can calculate either just the FFT or all three of the transformations. · For each transformation you can decide to look at the first n peaks in the signal. · And for each peak you can decide to take only the x value, or both the x and y values. In the example above, we have taken the x and y values of the first 5 peaks of each transform, so we have 270 columns (9*3*5*2) in total.

4、利用传统的scikit-learn库中的分类器进行分类

使用scikit-learn库中的随机森林分类器的结果如下。

from sklearn.ensemble import RandomForestClassifierfrom sklearn.metrics import classification_reportclf = 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))---Accuracy on training set is : 1.0---Accuracy on test set is : 0.9097387173396675---          1       0.96      0.98      0.97       496---          2       0.94      0.95      0.95       471---          3       0.94      0.90      0.92       420---          4       0.84      0.82      0.83       491---          5       0.86      0.96      0.90       532---          6       0.94      0.85      0.89       537------ avg / total      0.91      0.91      0.91      2947

但是，存在特征选择的问题：如果可以进一步地进行有效的特征选择，那么精度应该会有所提升。

It is understandable that some of the 270 features will be more informative than other ones. It could be that some transformations of some components do not have five peaks, or that the frequency value of the peaks is more informative than the amplitude value, or that the FFT is always more informative than the auto-correlation.The accuracy will increase even more if we actively select the features, transformations and components which are important for classification. Maybe we can even choose a different classifier or play around with its parameter values (hyperparameter optimization) to achieve a higher accuracy.

5、结语

随机信号分析为我们提供了一套对时间序列和信号进行分析、建模和分类的有力工具(FFT、PSD和自相关系数)。