持续创作,加速成长!这是我参与「掘金日新计划 · 10 月更文挑战」的第7天,点击查看活动详情
随机信号分析包括了信号的处理、分解和分析的科学领域。来自物理或者工程方面工程人员一定或多或少接触过信号,做过一些信号分析相关的工作,例如信号分析、信号建模和信号的分类。
术语表
- 频谱: 频谱是频率谱密度的简称,是频率的分布曲线。 复杂振荡分解为振幅不同和频率不同的谐振荡,这些谐振荡的幅值按频率排列的图形叫做频谱。
- 周期: 是周期信号瞬时幅值随时间重复变化的信号。周期性信号重复所需的时间称为周期 P
- 频率: 频率通常 f 是周期的倒数, 如果一个信号的周期为1秒,其频率为1赫兹,如果周期为10秒,其频率为0.1赫兹
- 赫兹:赫兹是国际单位制中频率的单位,表示每秒钟的周期性变动重复次数。指每秒发生周期波动的次数
- 傅里叶变换:
- 小波变换: 小波变换是新的一种变换分析方法,
- 奈奎斯特频率(Nyquist frequency) 是为防止信号混叠需要定义最小采样频率
随着计算机科学和统计学逐渐地渗透各个领域,专业人士也将这些技术应用到了信号分析。今天我们就来看一看如何使用将机器学习分类算法应用到随机信号分析中,从而准确地对时间序列和信号进行准确的分类和建模。
通过这篇文章,应该能够对一些信号处理技术有所认识,这些技术可用于从信号中提取特征,利用这些特征可以对信号进行,例如通过心电图来识别一个人。
信号相关基础知识
其实我们身边很多事物都是可以传达信息的信号,大家能够张口就来的信号,可能是手机信号、音频信号或者视频信号,这些信号都是一种形式来传递信息。除此之外图像也可以看作信号,图像可以看作 2 维空间上 RGB 亮度变化信息,还有例如雷达、声呐数据也是信号,还有就是医疗上脑电图、心电图和肌电图也是信号。
从此可见只要能够传递信息的东西,我们就可以将其视为信号。接下来我们在说一说信号和时间序列有什么区别
信号和时间序列
信号和时间序列到底是不是一个东西,当你被问到这个问题,第一感觉他们应该不是同一个东西吧,不过到底他们之间什么呢?可能不能立即给出答案。这里我来简单解释一下,时间序列是一个关于时间的函数,可以表示 y = y(t) 这样函数可以描述比特币价值或者股价随时间的变化,鱼群随着时间的变换也可以看作时间序列。而信号呢?信号这是一个更普遍的概念,在信号中 y 不一定是关于时间函数,可以是空间的函数,例如图像,例如 y = y(x,y) 关于 x 和 y 坐标的函数。
离散时间信号
通常信号处理会分为 2 个部分,若考虑离散时间,时域中的函数或信号,在离散时间点上的值是已知的,若考虑连续时间,在任意时间点上数值均是已知的。例如我们在 t=0.1,t=0.1435 或者 t=0.1435297 不断细分都可以得到对应值。
在自然界存在的大多数信号都是模拟信号的(连续的);想想我们人类体内的电信号,人类的语言,听到的任何其他声音,白天测量的光量,气压等等这些都是连续信号。
每当你想把这些模拟信号中的一个数字化,以便在计算机上对其进行分析和可视化,首先需要将信号变成了离散的。
我们可以看到,选择一个好的采样率是很重要的;如果采样率选得太低,离散信号将不再包含原始模拟信号的特征,定义信号的重要特征将丢失。
更具体地说,如果一个信号的采样率小于奈奎斯特率,就说是采样不足。奈奎斯特率是信号中最高频率的两倍。为了防止欠采样,通常选择比奈奎斯特率高得多的频率作为采样频率。
如果信号包含一个模式,在一个特定的时间段后重复出现,我们称之为周期性信号。
周期信号和傅里叶变换
如果信号模式在总是在一段时间重复出现,那么这样信号就称为周期信号
傅里叶分析是一个用于周期性分析的工具,或者说方法。如果一个信号包含周期性的成分,傅里叶分析可以用来将这些信号分解成周期性的成分。傅里叶分析
如果测量心跳,测量您心率为 60次/分钟,那么信号频率就是 1赫兹。每 2 秒钟移动一下手指,那么传到你手上信号就是 0.5 赫兹。如果将这个信号进行组合,然后再通过
在 python 中 scipy 库提供了关于傅里叶分析的方法,
由于我们的信号是以 100Hz 的频率 fsf_s 进行采样的,FFT 将返回频率为 f_s / 2 = 50Hz 的频谱。你的采样率越高,FFT能计算的最大频率就越高。
def 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_values
在上面的 get_fft_values
函数中,scipy.fftpack.fft
函数返回一个复值频率的向量。由于它们是复值,它们将包含一个实部和一个虚部。复数值的实部对应于信号的幅度,虚部对应于信号的相位。由于我们只对振幅的大小感兴趣,所以我们使用np.abs()来获取频谱的实部。
t_n = 10 #设置 10 秒
N = 1000 #采样的次数
T = t_n/N # 采样时间间隔
f_s = 1/T #采样频频率
x_value = np.linspace(0,t_n,N) #10 秒钟
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))]
for y_val in y_values:
plt.plot(y_val)
composite_y_value = np.sum(y_values,axis=0)
plt.plot(composite_y_value)
f_values,fft_values = get_fft_values(composite_y_value,T,N,f_s)
plt.plot(f_values,fft_values,linestyle='-',color='b')
plt.xlabel('Frenquency[HZ]',fontsize=16)
plt.ylabel("Amplitude",fontsize=16)
plt.show()
colors = ['k','b','b','b','b','b','b','b','b']
fig = plt.figure(figsize=(8,8))
ax = fig.add_subplot(111,projection='3d')
ax.set_xlabel("nTime[s]",fontsize=16)
ax.set_ylabel("nFrequency[Hz]",fontsize=16)
ax.set_zlabel("nAmplitude[s]",fontsize=16)
y_values_ = [composite_y_value] + list(reversed(y_values))
frequencies = [1,1,1.5,3,5,6.5]
print(frequencies)
for ii in range(0,len(y_values_)):
signal = y_values_[ii]
color = colors[ii]
length = signal.shape[0]
x = np.linspace(0,10,1000)
y = np.array([frequencies[ii]]*length)
z = signal
if ii == 0:
linewidth = 4
else:
linewidth = 2
ax.plot(list(x),list(y),list(z),linewidth=linewidth,color=color)
x = [10]*75
y = f_values[:75]
z = fft_values[:75]*3
ax.plot(list(x),list(y),list(z),linewidth=2,color='red')
plt.tight_layout()
plt.show()
功率谱密度
定义了信号或者时间序列的功率如何随频率分布。这里功率可能是实际物理上的功率,或者更经常便于表示抽象的信号被定义为信号数值的平方,也就是当信号的负载为 1 欧姆(ohm)的实际功率。由于平均值不为零的信号不是平方可积的,所以在这种情况下就没有傅里叶变换。 维纳-辛钦定理(Wiener–Khinchintheorem)提供了一个简单的替换方法,如果信号可以看作是平稳随机过程,那么功率谱密度就是信号自相关函数的傅里叶变换。
import numpy as np
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
plt.plot(f_values,psd_values,linestyle='-',color='blue')
plt.xlabel('Frequency [Hz]')
plt.ylabel("PSD[V**2]/Hz")
plt.show()