用Python使用快速傅里叶变换分析振动传感器采集的数据并绘制趋势图、分布图和频谱图

背景某振动传感器可以通过蓝牙将测量的设备振动信号传输到手机，现需要对采集到的数据进行分析，并绘制趋势图、数据分布图和频谱图。振动传感器的采样频率为12.8KHz（采样间隔为 1e6/12800=78.125微秒）,每秒钟最多可以将2048个（160ms的测量数据）数据传输到手机。采集获得的数据保存为文本文件，数据样式如下图。环境本文使用 python 3.9.6，在Windows 11环境下，用

bkzy

10046人浏览 · 2021-09-17 16:46:08

bkzy · 2021-09-17 16:46:08 发布

背景

某振动传感器可以通过蓝牙将测量的设备振动信号传输到手机，现需要对采集到的数据进行分析，并绘制趋势图、数据分布图和频谱图。
振动传感器的采样频率为12.8KHz（采样间隔为 1e6/12800=78.125微秒）,每秒钟最多可以将2048个（160ms的测量数据）数据传输到手机。采集获得的数据保存为文本文件，数据样式如下图。
在这里插入图片描述

环境

本文使用 python 3.9.6，在Windows 11环境下，用 jupyter notebook 进行分析。

读取数据文件

import numpy as np 
import pylab as pl #导入绘图模块
from matplotlib.pylab import mpl
mpl.rcParams['font.sans-serif'] = ['SimHei']   #显示中文
mpl.rcParams['axes.unicode_minus']=False       #显示负号

# 通过文件读取数据
with open("振动数据2.txt", "r") as f:  # 打开文件
    data = f.read()  # 读取文件
# 文件重点数据使用逗号分隔，此处需要把数据分隔提取
data = data.split(',')
# 提取到的数据是文本格式的，需要转换为浮点数的数据
x = [float(v) for v in data]
print("共有%d点数据"%(len(x)))
"""
此处输出为：
共有38912点数据
"""

绘制趋势图和数据分布直方图

sampling_rate = 12800 #取样频率(来自传感器说明书)
fft_size = 2048 #FFT处理的数据样本数

# 用linspace计算每个数据点的时间标，时间标的单位为ms,总时长为160ms(2048个数据点)：
# linspace在指定的间隔内返回均匀间隔的数字
# np.linspace(start, stop, num=50, endpoint=True, retstep=False, dtype=None)
tstemp = np.linspace(0,int(1e6/sampling_rate * fft_size),fft_size)/1e3

# 绘制图形
pl.figure(figsize=(14,8))
pl.title("原始数据(趋势图)")
pl.ylabel("振动加速度(m/s2)")
pl.xlabel("时间(ms)")
pl.plot(tstemp,x[:fft_size])

pl.figure(figsize=(14,8))
pl.title("原始数据(分布图)")
pl.ylabel("点数")
pl.xlabel("振动加速度(m/s2)")
pl.hist(x[:fft_size],bins=100)
pl.show()

绘制的图形如下：
在这里插入图片描述

执行傅里叶变换

在np.fft中，有fft和rfft可以选择。
如果样本数据的总数量为N，fft函数返回N个变换后的数据，但是这N个数据是左右对称的，即有效的实际数据只有N/2+1个。
而rfft直接只返回N/2+1个数据，可以不用处理直接使用。

# 进行快速傅里叶变换
xs = x[:fft_size]# 从波形数据中取样fft_size个点进行运算
# 利用np.fft.rfft()进行FFT计算
# rfft返回N/2+1个数据。而fft返回N个数据，但这N个数据是左右对称的，因此实际有效的是N/2+1个数据。
# 所以此处直接使用rfft
xfft = np.fft.rfft(xs) 
# 打印显示前5个傅里叶转换结果
# 可以看出傅里叶转换结果是复数
# 复数用y=a+bj的形式表示，a为实部，b为虚部
print(xfft[:5])

此处的输出数据为

[2442.331       +0.j         -216.7047702  -55.99042361j
 -394.94605398-121.14476753j    8.80615719 +90.98250801j
  -65.84171697 +10.04387751j]

绘制频谱图

取频率
绘制频谱图，首先一点试要知道横坐标上每点对应的频率值。
本项目中的采样频率sampling_rage=12800Hz，共有N=fft_size=2048个采样点。
将0~sampling_rage分为N份，然后将美分对应给各个数值，即使每个样本数据点所对应的频率。
由于经过rfft傅里叶变换后，频率和样本点折半了，因此计算中频率使用sampling_range/2,样本点数使用N/2+1。
复数取模运算
相当于始数取绝对值，但对于复数来说，没有绝对值这个概念。
复数取模运算是对实部和虚部的平方和开根号，即：
$\sqrt{a^2 + b^2}$

# 通过下面的np.linspace计算出返回值中每个下标对应的真正的频率：
# 最后除以1e3将频率单位由Hz转换为KHz
freqs = np.linspace(0,int(sampling_rate/2),int(fft_size/2+1))/1e3

# 对复数取绝对值(取模运算)，复数的模为：√a^2 + b^2
abs_xfft=abs(xfft)

pl.figure(figsize=(14,8))
pl.title("单边振幅谱(未归一化)")
pl.ylabel("振幅(未归一化)")
pl.xlabel("频率(KHz)")
pl.plot(freqs,abs_xfft)
pl.show()

绘制的图形为：
在这里插入图片描述
可以发现它的纵坐标非常大，看不出具体的物理意义。为了使纵坐标具有物理意义，工程应用中需要将取模后的傅里叶转换结果归一化，即将数值除以样本总数N.

数据归一化

# 直接对傅里叶转换后的结果进行归一化处理
# 归一化处理(可以在取模前归一化，也可以在取模后归一化)
abs_xfft_n = abs_xfft/fft_size

# 再次绘图
pl.figure(figsize=(14,8))
pl.title("单边振幅谱(归一化后)")
pl.ylabel("振动加速度(m/s2)")
pl.xlabel("频率(KHz)")
pl.plot(freqs,abs_xfft_n)
pl.show()

在这里插入图片描述
可以看到此时的图的纵坐标具备了一定的物理含义。但是最大值与最小值之间差别比较大，大值将小值压制的几乎只有一条线了，看不太清趋势。此时还可以对纵坐标的值取对数，缩小大值和小值之间的差别。但缺点是纵坐标的物理意义又不清楚了。

纵坐标取对数

# 对归一化取模后的数据再取以10为底的对数
# np.clip(a, a_min, a_max, out=None, **kwargs)
# np.clip将输入中a限制在min和max之间，将小于min的置为min，大于max的置为max
# 防止出现 log(0)导致无法计算的问题
xfft_log = np.log10(np.clip(abs_xfft_n,1e-20,1e1000))
pl.figure(figsize=(14,8))
pl.title("单边振幅谱(取Log10后)")
pl.ylabel("振动加速度(Log10(m/s2))")
pl.xlabel("频率(KHz)")
pl.plot(freqs,xfft_log)
pl.show()