马上注册,结交更多好友,享用更多功能^_^
您需要 登录 才可以下载或查看,没有账号?立即注册
x
本帖最后由 incrediblejack 于 2020-6-2 22:26 编辑
本例的主要功能:
将时域信号通过快速傅里叶变换FFT,转变为频域信号,并对频域信号提取特征,形成新的数据集。
其中,时域信号是很长的,傅里叶变换在原始时域信号上进行多次分段采样,之后的频域信号也是一段数据,而特征提取之后变成一个点,这样一段时域信号--对应一段频域信号--对应一个特征点,多段时域信号--对应多段频域信号--对应多个特征点。
最终形成一个集成了所有特征点的数据集excel文件。
但是目前有bug
错误为:
Traceback (most recent call last):
File "C:/Users/Administrator/PycharmProjects/feature202005/p_61_03.py", line 185, in <module>
ff1 = function_f1(xfp_value) #调用f_1函数 输入是加了绝对值的频谱幅值xfp_value
File "C:/Users/Administrator/PycharmProjects/feature202005/p_61_03.py", line 47, in function_f1
sum_mid_f1_v1 = input[i]
IndexError: index 125 is out of bounds for axis 0 with size 125
源代码已附上。
各位大神帮忙看看咯!
谢谢!
另外,怎么添加附件呀?我没找到添加附件的选项,我的框框中只有表情和图片,没有“文件”那个按钮啊,所以我附上个链接吧,里面是数据集。
链接: https://pan.baidu.com/s/1U4gADju6QtH2ZRFDLEb3-g 提取码: fgcu
#这个版本是在p_61_02基础上改的。
#20200612 19:16
#时域和频域分别绘图
import numpy as np
import matplotlib.pyplot as plt
import os
os.environ['TF_CPP_MIN_LOG_LEVEL'] = '2'
import xlrd
import math
import xlwt
workbook = xlwt.Workbook(encoding='utf-8') #写入数据是用的临时变量
booksheet = workbook.add_sheet('Sheet 1', cell_overwrite_ok=True)
train1 = xlrd.open_workbook('shuangtu_13_14.xlsx') #读取原始数据文件
table1 = train1.sheets()[0] #提第一页sheet
len_train = 1000 #设置取用数据长度(可更改)
# p_a = [0 for i in range(len_train)] #中间变量,列表p_a
# for i in range(len_train): #1000次
# p_a[i] = table1.row_values(i)[1] #将第1列的内容,逐行读取,存入p_a列表
#原始数据集装入a
a = [0 for i in range(len_train)] #中间变量,列表a
for i in range(len_train): #1000次
a[i]=table1.row_values(i)[1] #将第1列的内容,逐行读取,存入a列表
fft_size = 250 #快速傅里叶变换长度
start_i = 1 #设置初始提取位置
#傅里叶变换的方法,封装成函数:
def Fast_Fourier_Trans(x):
sampling_rate = 500 #采样频率
t = np.linspace(0, 1, sampling_rate) # 截取一段时间,截取是任意的,这里截取了0~1秒的一段时间。
xs = x[start_i:start_i + fft_size] #进行原始数据分割和提取,将用于进行傅里叶变换
xf = 2*np.fft.rfft(xs) / fft_size # 傅里叶变换公式,得到频谱幅值
freqs = np.linspace(0, sampling_rate // 2, fft_size // 2 + 1) #横轴: 表示频率, 从0开始,到采样频率//2,一共显示有包括首尾的fft_size/2 +1个数
xfp = np.clip(np.abs(xf), 1e-20, 1e100) #为了避免频谱幅值出现负数,添加的取绝对值的步骤。
t_heng = t[start_i:start_i + fft_size] # 为了后面显示用
return (freqs, xfp, t_heng, xs, t, xf) #多个返回值。
#特征函数,求f_1 对应公式5-12
def function_f1(input):
sum_mid_f1_v1 = 0
sum_mid_f1_v2 = 0
for i in range(len(freqs_value)):
sum_mid_f1_v1 = input[i]
sum_mid_f1_v2 += sum_mid_f1_v1
f_1 = sum_mid_f1_v2 / len(freqs_value)
# print('平均能量,即f1:')
# print(f_1)
# print(' ')
return f_1
#特征函数,求f_2 对应公式5-13
def function_f2(input):
sum_mid_f2_v1 = 0
sum_mid_f2_v2 = 0
for i in range(len(freqs_value)):
sum_mid_f2_v1 = (input[i] - function_f1(input)) ** 2 # 式子5-13分子括号里的
sum_mid_f2_v2 += sum_mid_f2_v1 # 式子5-13分子。
f_2 = sum_mid_f2_v2 / (len(freqs_value) - 1)
# print('f2:')
# print(f_2)
# print(' ')
return f_2
#特征函数,求f_3 对应公式5-14
def function_f3(input):
sum_mid_f3_v1 = 0
sum_mid_f3_v2 = 0
for i in range(len(freqs_value)):
sum_mid_f3_v1 = (input[i] - function_f1(input)) ** 3
sum_mid_f3_v2 += sum_mid_f3_v1
f_3 = sum_mid_f3_v2 / (len(freqs_value) * ((function_f2(input)) ** 1.5))
# print('f3:')
# print(f_3)
# print(' ')
return f_3
#特征函数,求f_4 对应公式5-15
def function_f4(input):
sum_mid_f4_v1 = 0
sum_mid_f4_v2 = 0
for i in range(len(freqs_value)):
sum_mid_f4_v1 = (input[i] - function_f1(input)) ** 4
sum_mid_f4_v2 += sum_mid_f4_v1
f_4 = sum_mid_f4_v2 / (len(freqs_value) * ((function_f2(input)) ** 2))
# print('f4:')
# print(f_4)
# print(' ')
return f_4
#特征函数,求f_5 对应公式5-16
def function_f5(input):
sum_mid_f5_v1 = 0
sum_mid_f5_v2 = 0
sum_vid_f5_v3 = 0
sum_mid_f5_v4 = 0
for i in range(len(freqs_value)):
sum_mid_f5_v1 = (freqs_value[i]) * (input[i])
sum_mid_f5_v2 += sum_mid_f5_v1
sum_mid_f5_v3 = input[i]
sum_mid_f5_v4 += sum_mid_f5_v3
f_5 = sum_mid_f5_v2/ sum_mid_f5_v4
# print('f5:')
# print(f_5)
# print(' ')
return f_5
#特征函数,求f_6 对应公式5-17
def function_f6(input):
sum_mid_f6_v1 = 0 # fk^2*X(k)
sum_mid_f6_v2 = 0 # ∑fk^2*X(k)
sum_mid_f6_v3 = 0 # (∑fk^2*X(k))^0.5
sum_mid_f6_v4 = 0 # 分母
for i in range(len(freqs_value)):
sum_mid_f6_v1 = ((freqs_value[i]) ** 2) * (input[i])
sum_mid_f6_v2 += sum_mid_f6_v1
sum_mid_f6_v3 = (sum_mid_f6_v2) ** 0.5 # 这个好像算不出来。
sum_mid_f6_v4 += input[i]
f_6 = sum_mid_f6_v3 / sum_mid_f6_v4
# print('f6:')
# print(f_6)
# print(' ')
return f_6
#特征函数,求f_7 对应公式5-18
def function_f7(input):
sum_mid_f7_v1 = 0 # fk^4*X(k)
sum_mid_f7_v2 = 0 # ∑fk^4*X(k)
sum_mid_f7_v3 = 0 # (∑fk^4*X(k))^0.5
sum_mid_f7_v4 = 0 # fk^2*X(k)
sum_mid_f7_v5 = 0 # ∑fk^2*X(k)
for i in range(len(freqs_value)):
sum_mid_f7_v1 = ((freqs_value[i]) ** 4) * (input[i])
sum_mid_f7_v2 += sum_mid_f7_v1
sum_mid_f7_v3 = (sum_mid_f7_v2) ** 0.5
sum_mid_f7_v4 = ((freqs_value[i]) ** 2) * (input[i])
sum_mid_f7_v5 += sum_mid_f7_v4
f_7 = sum_mid_f7_v3 / sum_mid_f7_v5
# print('f7:')
# print(f_7)
# print(' ')
return f_7
#特征函数,求f_8 对应公式5-19
def function_f8(input):
sum_mid_f8_v1 = 0 # fk^2*X(k)
sum_mid_f8_v2 = 0 # ∑fk^2*X(k)
sum_mid_f8_v3 = 0 # fk^4*X(k)
sum_mid_f8_v4 = 0 # ∑fk^4*X(k)
sum_mid_f8_v5 = 0 # X(k)
sum_mid_f8_v6 = 0 # ∑X(k)
sum_mid_f8_v7 = 0 # (∑fk^4*X(k))*(∑X(k))
for i in range(len(freqs_value)):
sum_mid_f8_v1 = ((freqs_value[i]) ** 2) * (input[i])
sum_mid_f8_v2 += sum_mid_f8_v1
sum_mid_f8_v3 = ((freqs_value[i]) ** 4) * (input[i])
sum_mid_f8_v4 += sum_mid_f8_v3
sum_mid_f8_v5 = input[i]
sum_mid_f8_v6 += sum_mid_f8_v5
sum_mid_f8_v7 = (sum_mid_f8_v4) * (sum_mid_f8_v6)
f_8 = (sum_mid_f8_v7) ** 0.5
# print('f8:')
# print(f_8)
# print(' ')
return f_8
#主流程
for j in range(len_train // fft_size): #这个变量是为了控制写入数据的行数 j
if j < 2:
for er in range(len_train // fft_size): #这个变量是为了控制运行傅立叶变换Fft函数的次数 er
Fft = Fast_Fourier_Trans(a) #调用快速傅里叶变换函数Fft,输入是提取了len_train长度的原始a列表数据
freqs_value = Fft[0] #Fft返回值,这个是干啥的?
xfp_value = Fft[1] #Fft返回值,返回加了绝对值的频谱幅值
t_heng_value = Fft[2] #Fft返回值,返回当前时域(Fft变换之前的)横坐标值
xs_value = Fft[3] #Fft返回值,返回当前时域(Fft变换之前的)纵坐标值
t_value = Fft[4] #Fft返回值,返回采样间隔么?
xf_value = Fft[5] #Fft返回值,返回未加绝对值的频谱幅值(可能会有负的,所以不用这个)
ff1 = function_f1(xfp_value) #调用f_1函数 输入是加了绝对值的频谱幅值xfp_value
print("ff1:")
print(ff1) #测试用
print(' ')
ff2 = function_f2(xfp_value) #调用f_2函数。。。下面道理类同
print("ff2:")
print(ff2)
print(' ')
ff3 = function_f3(xfp_value)
ff4 = function_f4(xfp_value)
ff5 = function_f5(xfp_value)
ff6 = function_f6(xfp_value)
ff7 = function_f7(xfp_value)
ff8 = function_f8(xfp_value)
rowdata = [ff1, ff2, ff3, ff4, ff5, ff6, ff7, ff8] #形成一行数值
print('rowdata:')
print(rowdata) #测试用
print(' ')
for i in range(len(rowdata)): #这个变量i是为了遍历上面的rowdata,为了将来写入excel表格
booksheet.write(j,i, rowdata[i]) #读取rowdata的第i个,写入excel的第j列,第i行
i += 1
er += 1
start_i += fft_size #执行完一轮Fft之后变量start_i增加一个采样周期fft_size之后重新计算
j += 1
else:
break
plt.figure(figsize=(8,8)) #生成画布
plt.subplot(211) #生成子图 #这个图是
plt.plot(t_heng_value, xs_value) #时域图
plt.xlabel(u'时间(秒)', fontproperties='FangSong')
plt.title(u'原始幅值',fontproperties='FangSong')
plt.subplot(212)
plt.plot(freqs_value, xfp_value) ##频域图
plt.xlabel(u'频率(HZ)', fontproperties='FangSong')
plt.ylabel(u'幅值',fontproperties='FangSong')
plt.subplots_adjust(hspace=0.4)
plt.show()
# print('xfp_value')
# print(xfp_value)
# print(' ')
# print('xfp_value 傅立叶变换后的点纵坐标的个数:')
# print(len(xfp_value))
# print(' ')
# print('xs_value 傅里叶变换后的点横坐标的个数:')
# print(len(xs_value))
# print(' ')
workbook.save('p_61_02.xls') #保存
|