平时工作中测试NVH，一般使用PAK统来进行数据的采集和分析，但由于软件license限制原因，在不同的电脑上工作会很不方便。
这里介绍下如何用python来处理振动数据。数据来源于一次固有频率测试，信号接收器为麦克风，数据采集设备为MKII，用PAK导出了原始的时间振动采集数据，采样频率为102400赫兹。

首先是导入模块与文件：

'''Fourtran'''
'''Author:jAEgerrr'''
'''2020-02-28'''
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
%matplotlib notebook

df=pd.read_table('./throughput.txt',header=None,index_col=0)
df

通过read文件发现，文件有两列：

其中第0列为时间数据，作为index，第1列为振动数据（声压）。
我们首先对原始数据直接进行plot，看下时域上的振动信号:

plt.figure()
plt.plot(df,color='r')
plt.xlabel('Time[s]')
plt.ylabel('Sound Pressure[Pa]')
plt.title('Time Resopnse')

得到了一个声压随时间变化的图像，接下来我们想看一下这个信号在频率上有哪些特征，所以要进行傅里叶变换。
代码如下：

x_time=df.index
y_amp=df[1].values
rate=102400

'''FFT operation, shift zero frequency to the middle of the array'''
y_f=np.fft.fftshift(np.fft.fft(y_amp))
x_f=np.fft.fftshift(np.fft.fftfreq(x_time.size,d=1/rate))

''' APS spectrum, and switch Pa to dB '''

y_f_norm=np.sqrt(np.abs(y_f**2))
y_log=20*np.log10(y_f_norm)

'''Remove negative frequency'''

x_f_h=x_f[range(x_f.size//2+1,x_f.size)]
y_log_h=y_log[range(y_log.size//2+1,y_log.size)]  

'''Figure Operation'''
plt.figure()
plt.plot(x_f_h,y_log_h,color='k')
plt.xlabel('Frequency[Hz]')
plt.ylabel('Amplitude[dB]')
plt.title('Frequency Response')   

其中涉及傅里叶变换的地方：
np.fft.fft(y_amp)是直接把时域y轴数据（时间vs声压）变换成频域的y轴数据（频率vs声压）
np.fft.fftfreq(x_time.size,d=1/rate) 把时域的x轴数据（时间）变幻成了频域的x轴数据（频率）
np.fft.fftshift是把傅里叶变换后的数据以0Hz作为中心点从小到大排列，因为傅里叶变换后的数据有正频率也有负频率，需要排序一下
y_f_norm=np.sqrt(np.abs(y_f**2))
y_log=20*np.log10(y_f_norm）
是求y的自功率谱，然后把声压转换成分贝形式。
变换后的频谱图显示如下：

可以看到在16176Hz附近有一个明显的峰值，这个就是我们感兴趣的固有频率特征。
大功告成~~