信号处理技术在机器学习中的应用（二）

这是一个学习总结。

原文链接：

Machine Learning with Signal Processing Techniques – ML Fundamentals

2.时域和频域之间的转换

¢ 傅里叶变换把信号中的不同周期的组成部分分离出来

¢ 傅里叶变换可以告诉我们一个信号的每个具有周期的组成部分的频率

¢ 例如，如果你测量自己的心率，假设是 60 次 / 分钟，那么心跳的的周期是 1s ，则频率是 1Hz 。与此同时，你每两秒动一下你的手指头，你的手传出来的信号的频率则是 0.5Hz 。在你手臂上的电极测量了的信号是这两种信号的组合。傅里叶变换可以把这个合成的信号分离。你会看到这个合成的信号经过傅里叶变换后再频谱上有两个峰，一个峰在 0.5Hz ，一个峰在 1Hz 。

¢ 所有复杂的信号都可以通过傅里叶变换分解成一堆更加简单的信号

¢ 这些更加简单的信号就是三角函数（ sin 波和 cos 波）

¢ 傅里叶变换把信号从时域转换到频域

¢ 反傅里叶变换把信号从频域转换到时域

5 个 sin 波 （蓝色），振幅分别是 4 ， 6 ， 8 ， 10 和 14 ，频率分别是 6.5 ， 5 ， 3 ， 1.5 以及 1Hz
这 5 个信号合成一个 信号 （黑色）
傅里叶变换可以把 合成信号 转换到 频域 （红色）

 上图的绘制方法如下代码所示

from mpl_toolkits.mplot3d import Axes3D

t_n = 10
N = 1000
T = t_n / N   #合成波的周期
f_s = 1/T   #合成波的频率

x_value = np.linspace(0,t_n,N)#产生N个数，x_value是一个向量，向量属于numpy的数组，向量不能被 
                              #转置,矩阵也属于numpy数组，但是矩阵可以被转置

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))]  
#y_values是一个含有5个元素的列表，每个元素是一个numpy数组（即一个含有N个点的向量）
#注意：y=Asin(wx+r)+u
#r和u是在x轴和y轴上的位移
#A是振幅，w决定周期，最小周期T=2pi/|w|，因此编程这里|w|=2pi/T，w=2pi*Frequency

composite_y_value = np.sum(y_values, axis=0) #合成波，是一个含有N个数的向量

f_values, fft_values = get_fft_values(composite_y_value, T, N, f_s)

colors = ['k', 'b', 'b', 'b', 'b', 'b', 'b', 'b', 'b']

fig = plt.figure(figsize=(8,8))
ax = fig.add_subplot(111, projection='3d')#创建3d图形
ax.set_xlabel("\nTime [s]", fontsize=16)
ax.set_ylabel("\nFrequency [Hz]", fontsize=16)
ax.set_zlabel("\nAmplitude", fontsize=16)

y_values_ = [composite_y_value] + list(reversed(y_values))
frequencies = [1, 1.5, 3, 5, 6.5]

for ii in range(0,len(y_values_)):
    signal = y_values_[ii] #5个sin波
    color = colors[ii]  #sin波的颜色
    length = signal.shape[0] #波的长度，即向量的长度1000
    x=np.linspace(0,10,1000)  #x轴
    y=np.array([frequencies[ii]]*length) #y轴
    z=signal #z轴

    if ii == 0:
        linewidth = 4
    else:
        linewidth = 2
    ax.plot(list(x), list(y), zs=list(z), linewidth=linewidth, color=color)#时域

    x=[10]*75
    y=f_values[:75]
    z = fft_values[:75]*3
    ax.plot(list(x), list(y), zs=list(z), linewidth=2, color='red')
    
    plt.tight_layout() #tight_layout会自动调整子图参数,使之填充整个图像区域。
plt.show()

Python中的快速傅里叶变换

¢ 快速傅里叶变换（ FFT ）能够有效率地计算离散傅里叶变换。 FFT 是一个计算傅里叶变换的业界标准。 很多库、包都 包含 FFT 。

¢ python 中用来计算信号的 FFT 的包是 SciPy

from scipy.fftpack import fft

def get_fft_values(y_values, T, N, f_s):
    f_values = np.linspace(0.0, 1.0/(2.0*T), N//2)  #频率，最高为50Hz
    fft_values_ = fft(y_values)
    fft_values = 2.0/N * np.abs(fft_values_[0:N//2])
    return f_values, fft_values

t_n = 10
N = 1000
T = t_n / N
f_s = 1/T

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

 

合成波经过FFT后，从时域变换到频域

因为我们的合成波的频率fs是100Hz，

FFT返回频谱的范围的上限是fs/2,即50Hz。

sciPy包的fft会返回一个包含复数的向量，这个复数的实部是信号的振幅，虚部是信号的相位。

因为我们只对振幅感兴趣，使用np.abs取出复数的实部。

信号有N个点，则FFT后，仍然返回N个点。

其中前一半的点N/2含有有用的值，这些值对应频谱的范围从0Hz到Nyquist rate即fs/2。

注意：在y=Asin(ωx+φ)中，A称为振幅；ωx+φ称为相位；x=0时的相位（ωx+φ=0+φ=φ）称为初相。

Python中的功率频谱密度PSD

PSD即Power Spectral Density（功率频谱密度）

¢ 跟 FFT 类似，用来描述信号的频谱

¢ 还 会考虑信号每个频率的功率分布

¢ 一般来说，频率的峰的位置跟 FFT 的一样，但是峰的高度和宽度跟 FFT 的不一样

¢ 每个峰下面的面积表示频率的分布

from scipy.signal import welch

def get_psd_values(y_values, T, N, f_s):
    f_values, psd_values = welch(y_values, fs=f_s)
    return f_values, psd_values


t_n = 10
N = 1000
T = t_n / N
f_s = 1/T

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

 用python计算自相关

¢ 如果一个信号里面有一个经过 t 时间重复出现的模式，那么这个信号与延迟 t 时间的版本的信号具有高度自相关性

¢ PSD 可以通过 FFT 自相关函数得到，自相关函数可以通过把 PSD 函数 IFFT 得到

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_values

t_n = 10
N = 1000
T = t_n / N
f_s = 1/T

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

 小波变换

¢ 小波变换和傅里叶变换一样，把信号从时域转换到频域

¢ 傅里叶变换只有频域的高度解析，没有时域的解析。即我们可以知道频率在哪个位置振荡，但是不知道这个振荡在什么时候发生。

¢ 小波变换无论在频域还是时域都有高度的解析。

¢ 小波变换适合用来分析具有动态频谱的信号。即频谱随着时间改变会发生 变化