使用Jupyter Notebook进行波形分析
首先拿到了振动传感器的数据
解释数据:
我们知道机械波分为横波
、竖波
、以及复杂波
,现实世界中全部为复杂波——上下左右都有振动。
传感器每秒记录10240条数据(**每分钟仅记录前30秒的数据**
),对应X轴
、Y轴
、以及Z轴
的数据。我们这次仅分析一天
Z轴
的数据。
由于一秒10240条数据一天的数据量就是
10240*30s*60min*24h = 442368000条记录
,这个数据量挺大的,如果将一天的数据都画出波形图,会导致没有任何细节,整个图都会缩在一起。所以我这里先做了每分钟的均方根
处理,然后每分钟去画图能够清晰的展示出每天哪个时间段有较强的振动。
--------------------------------接下来看代码以及效果----------------------------
#以下是导入的包、以及用到的方法
import numpy as np
from scipy.fftpack import fft,ifft
import matplotlib.pyplot as plt
import os
import math
#计算均方根
def get_rms(records):
"""
均方根值 反映的是有效值而不是平均值
"""
return math.sqrt(sum([x ** 2 for x in records]) / len(records))
#读取文件
def readname():
filePath = '/root/Rone/waveData/data'
name = os.listdir(filePath)
return name
接着将一天的数据上传到服务器指定的文件夹 ——计算每分钟的均方根,以及将结果保存(这里我打印了出来,个人感觉比保存到文件中更方便去处理)
#求每分钟的RMS
names = readname()
arr = []
for i in names:
tempData = np.loadtxt("/root/Rone/waveData/data/"+i)
rmsData = get_rms(tempData)
arr.append(rmsData)
print(i)
file=open('/root/Rone/waveData/resultRMS.txt','w')
file.write('\n'.join(str(arr)))
file.close()
print("=============")
for i in arr:
print(i)
下图是部分均方根的执行结果
使用matlab库画出一天的均方根图
dataRMS = np.loadtxt("/root/Rone/waveData/resultRMS.txt")
x = np.linspace(0,1313,1313)
plt.subplot(4,1,1)
fig=plt.gcf()
fig.set_size_inches(128.5, 100.5)
plt.plot(x,dataRMS[0:1313])
plt.title('24h Data Analysis')
plt.xlabel('t')
plt.ylabel('x')
x2 = np.linspace(500,1313,813)
#print(x2)
print(dataRMS.size)
fig=plt.gcf()
fig.set_size_inches(20.5, 10.5)
结果图(16点——第二天14点的数据)我们可以看出峰值高的是下午5点左右 以及第二天8-10点左右。(使用均方根后出现一个问题——没有了负值,这是因为均方根是经过平方后再开方导致的)。
均方根只会把差距放大,不能代表波形的振幅
。
接下来我们就去找异常值。
#查询16点到16.30点的数据 total
names = readname()
arr = np.zeros(1)
for i in names:
tempData = np.loadtxt("/root/Rone/waveData/data/"+i)
#
#print(tempData)
arr = np.concatenate((arr,tempData),axis=0)
if i =='20190909_165228.txt':#
break;
print(i)
print(len(arr))
lenthData = len(arr)
np.savetxt(r"/root/Rone/waveData/16_17_hour.txt",arr)
x = np.linspace(0,lenthData,lenthData)
#print(x)
plt.subplot(4,1,1)
plt.plot(x,arr[0:])
fig=plt.gcf()
fig.set_size_inches(158.5, 130.5)
#查询16.30点到17.00点的数据 total
names = readname()
arr = np.zeros(1)
rootPath = "/root/Rone/waveData/data/"
# 52 58.txt
for i in (53,59):
tempData = np.loadtxt(rootPath+"20190909_16"+str(i)+"58.txt")
#
#print(tempData)
arr = np.concatenate((arr,tempData),axis=0)
if i =='20190909_170058.txt':#
break;
print(i)
print(len(arr))
lenthData = len(arr)
#np.savetxt(r"/root/Rone/waveData/16_30_17_00_hour.txt",arr)
x = np.linspace(0,lenthData,lenthData)
plt.subplot(4,1,1)
plt.plot(x,arr[0:])
fig=plt.gcf()
fig.set_size_inches(158.5, 130.5)
#查询17.00点到17.30点的数据 total
names = readname()
arr = np.zeros(1)
rootPath = "/root/Rone/waveData/data/"
# 52 58.txt
for i in range(0,9):
tempData = np.loadtxt(rootPath+"20190909_170"+str(i)+"58.txt")
arr = np.concatenate((arr,tempData),axis=0)
if i =='20190909_170958.txt':#
break;
print(i)
for i in range(10,30):
#20190909_172958
tempData = np.loadtxt(rootPath+"20190909_17"+str(i)+"58.txt")
arr = np.concatenate((arr,tempData),axis=0)
if i =='20190909_170058.txt':#
break;
print(i)
print(len(arr))
lenthData = len(arr)
#np.savetxt(r"/root/Rone/waveData/16_30_17_00_hour.txt",arr)
x = np.linspace(0,lenthData,lenthData)
#print(x)
plt.subplot(4,1,1)
plt.plot(x,arr[0:])
fig=plt.gcf()
fig.set_size_inches(158.5, 130.5)
上图有断点数据 很大的可能是因为传感器只有每分钟前30秒的数据。
接着查询
#查询21.40点到22.50点的数据 total
names = readname()
arr = np.zeros(1)
rootPath = "/root/Rone/waveData/data/"
for i in range(40,41):
tempData = np.loadtxt(rootPath+"20190909_21"+str(i)+"58.txt")
arr = np.concatenate((arr,tempData),axis=0)
print(i)
#lenthData = len(arr)
lenthData = 10240
x = np.linspace(0,1,10240)
#原始波形
plt.subplot(4,1,1)
plt.plot(x,arr[0:10240])
#加窗
w = np.hanning(10240)
plt.subplot(4,1,2)
plt.plot(x,w)#画出 hanning窗
z = w*arr[0:10240]
plt.subplot(4,1,3)
plt.plot(x,z)#加hanning窗后的波形
fig=plt.gcf()
fig.set_size_inches(30.5, 16.5)
接着进傅里叶转换
yf1=abs(fft(z[0:])) #归一化处理
yf2 = yf1[range(5120)] #由于对称性,只取一半区间
xf = np.arange(5120)
plt.subplot(222)
plt.plot(xf[0:],yf2[0:])
plt.title('FFT',fontsize=12)
fig=plt.gcf()
fig.set_size_inches(158.5, 130.5)
本篇到此结束,下一篇我分享一下低通滤波后的结果+傅里叶展开后的特征图形。