背景及摘要很多初学者刚接触故障诊断可能觉得很简单，套用深度学习模型进行训练，分类准确率很容易就达到99%。但是在工作或者在解决实际问题时，深度学习方法几乎是用不上（至少目前是这样）。 最近在读时献江老师等人的《机械故障诊断及典型案例解析》这本书时，发现其内容真的很丰富，故障诊断方法讲解的详细全面，且有实际丰富的案例，更加容易消化理解各方法特点及其应用场景。 这期给大家分享时域信号是如何转换到频域中的，在实际应用中又是如何将离散信号转换成频谱的。通过阅读本文，你将会学到傅里叶级数、傅里叶变换、采样定理等知识，还会了解到什么是频率混叠、频率泄露、栅栏效应。 参考资料：书籍：机械故障诊断及典型案例解析（第2版，时献江） 目录 1. 傅里叶级数 1.1 傅里叶级数原理及典型信号频谱特点2. 典型信号的傅里叶变换 2.1 单位脉冲信号（ 函数） 2.1 周期单位脉冲信号（梳妆 函数）3. 信号的采样 3.1 连续信号的采样 3.2 采样定理 3.3 采样点数和频率分辨率4. 离散傅里叶变换（DFT） 4.1 DFT的理论公式及计算过程5. 快速傅里叶变换（FFT） 5.1 FFT方法提出背景 5.2 FFT计算案例摘要随机信号的时域分析只能提供有限的时域特征信息，若想进行信号的精密分析与诊断，需要更进一步的频谱分析。因为故障发生时往往会引起信号频率成分的变化，而且很多故障特征频率也是可以计算的，通过监测该频率的幅值变换规律，就可以监控故障的发展过程。频谱分析的理论基础是傅里叶变换，傅里叶变换包括傅里叶级数和傅里叶积分。本章主要介绍傅里叶级数及典型信号的傅里叶变换、采样定理以及DFT的过程。 通过阅读本文将了解到（或带着以下疑问去阅读） 频域是什么，频域包含哪些信息？时域是怎么转换到频域的，时域转换到频域的原理是什么？连续信号转换到离散信号有什么要注意的吗？ FT、DFT、FFT之间是什么关系？时间分辨率和频率分辨率之间的关系？直流分量是什么，会有什么影响？ 1、傅里叶级数1.1 傅里叶级数原理及典型信号频谱特点(1) 基本原理给定一个周期函数 ，在一定条件下，可以根据如下公式展成傅里叶级数。 (1)式中 ——基频，Hz， ——周期，s。系数 和 用下面的公式进行计算 ，即： (2) (3) (4)式中， 是 的均值，称为直流分量， ， 为交流分量，也称谐波分量。在数学上，傅里叶级数的概念是任何一个周期函数都可以用无限个三角函数去拟合。在信号处理领域的物理意义为：对任何一个周期信号都可以分解成直流成分和无限个不能再分解的简谐信号（谐波分量）的叠加，或称任意周期信号中包括的频率成分如式(2)和式(3)所示。工程信号处理的一个主要任务就是分析和处理复杂信号中所包含的频率成分，这可以由傅里叶变换来完成。利用和差化积公式，我们可以得到更简洁的傅里叶级数表达式： (5)于是，幅值 和相位 ，可分别表示为： (6) (7)式中， 也可以采用复数描述傅里叶级数，根据欧拉公式有： (8) (9)代入式 (1) 式变为： 即： (10)式中， 。此时出现了负频率，这是由于用复数表示 ， 函数的结果。下面求系数 ： 可见 , 是一个复数量。由于 本身可以表示信号的幅值和相位，所以 称为有限区间 上信号 的离散频谱。 其中， 称为幅值谱，为复数谱的幅值，为实谱幅值 的一半；， 称为相位谱。(2) 周期信号傅里叶级数谱的特点根据式 (5) 可以绘出傅里叶级数的幅值谱图如图1所示。可见，周期信号的频谱是离散的，每条谱线只出现在基波频率的整倍数上，不存在非整倍数的频率分量，随着谐波次数的增加，谐波幅值是逐渐下降的。根据这些特征，很容易在频域内判别信号是周期还是非周期的。图1 傅里叶级数的幅值谱图例如，图2(c)中的两个周期信号的频率分别为2Hz和3Hz, 其合成显然为周期信号，且基频为1Hz 。图2(d)为其计算频谱图，并没有图1的特征，其实，此时如果认为基频幅值为零，其仍然具有周期函数的主要频谱特征。另外，也可对图3(c)中的合成信号进行频谱分析，结果如图3(d)所示，显然图中两个信号频率不成比例关系，虽然其时域波形很像周期信号，但是其谱特征为非周期信号（实际为准周期信号）。图2 两个正弦信号的叠加(有公共周期)图3 两个正弦信号的叠加(无公共周期)(3) 典型信号的傅里叶级数谱傅里叶级数的典型例子是方波信号[见图4(a)]，方波的傅里叶级数表达式为： 图4 方波及其频谱方波的频谱特点是只有奇次谐波。根据方波信号的频谱图4(b) ，也很容易判断其为周期信号。2、典型信号的傅里叶变换2.1 单位脉冲信号 ( 函数） (1) 函数的定义 (11)如图5(b)所示， 函数可以看作是图5(a)所示的矩形函数当 的极限。 图5 单位脉冲函数(2) 函数的性质① 函数的筛选性质a. 任何一函数 与 相乘的积分值等于此函数在零点的函数值 。 b. 任一函数 与具有时移 的单位脉冲函数 乘积的积分值是在该时移点上此函数的函数值 ： 根据这些性质可以看到函数 与处于时间轴上某点的单位脉冲函数相乘后的积分都等于该点的函数值，而其他所有点皆为零。这样，就等于通过这一处理将任意点的函数值筛选出来。这一特性后面被用来描述信号的采样过程。 ② 函数的卷积性质a.任一函数 与 的卷积仍是此函数本身: (12)b.同理，任一函数 与具有时移 的单位脉冲函数 的卷积是时移后的该函数 ： (13)也可以采用图6所示的图解方法描述，一个在坐标原点对称的矩形函数 与具有时移的两个单位脉冲 和 作卷积，其结果是将 移至在这两个单位脉冲函数所在位置上，即一函数与时间轴上任一点的单位脉冲函数卷积，其结果是将该函数原封不动地搬移到单位脉冲函数所在的时间轴位置。 图6 任一函数与 函数的卷积(3) 函数的频谱 函数时频域波形及对应频谱如图7所示，可知其在整个频域上都有分布。 图7 函数的频谱(4) 函数的物理意义 函数也叫冲击函数，是用以把一些抽象的不连续的物理量表示成形式上连续且能进行各种数学运算的广义函数。例如，它可以把集中载荷表示为分布载荷（分布面积为零，载荷集度为无穷大），把理想的碰撞冲量表示成一般冲量（作用的时间为零，作用力为无穷大 ），可以作为一种理想的脉冲激励函数。例如，采用敲击实验法测量叶片的共振频率，就是利用脉冲激励原理实现的。由锤敲击产生的冲击脉冲相当一个 函数，其频谱在理论上是一根直线，也就是在所有的频段具有相同的能量，称为白噪声。在与叶片的固有频率交叉处，会产生共振，引起系统的自由衰减振动，检测这个振动频率就是被测叶片的固有频率。 反之，直流信号的频谱是 函数，所以当信号中含有直流成分时，在频谱中的0频率附近就会有较高幅值的谱峰 ( 函数）出现，这会严重影响其他计算结果的显示比例，而且处于0点附近，不易被察觉，易产生误判，如图8(a) 所示。所以，如前所述，信号处理前一般都要进行均值化处理。零均值化处理后的频谱如图8(b) 所示，可见，简谐信号的幅值已变为最大值，说明 函数的影响被剔除。图8 含有均值和不含均值的简谐信号的频谱 2.2 周期单位脉冲信号 (梳状 函数）(1) 定义周期单位脉冲序列 如图9(a)所示，其解析表达式为: (14)式中， 为周期脉冲序列 的周期。图9 周期单位脉冲序列及其频谱(2) 梳状 函数的频谱频谱图如9(b)所示，周期单位脉冲序列的频谱仍是一个周期单位脉冲序列，其幅值为时域周期 的倒数， 即 ；频谱的周期也是时域周期的倒数。梳状函数主要用于从数学上描述信号的采样过程。3、信号的采样 在信号处理领域，计算机只能处理离散的时间序列函数，因此计算机需要把连续变化的信号变成离散信号后再进行相关的处理与运算。3.1 连续信号的采样一个连续时间信号 的采样过程如图10所示，图10(b)中的采样开关称为采样器，在采样器作用下的实际采样信号用 表示。信号处理系统的采样一般多为定时等间隔周期采样，采样间隔为 ，此时，当 时刻，采样开关闭合，采样器的输出恒等于 采样开关闭合时间为 ，到达 时刻，采样开关打开，采样器的输出为零。这样，在采样器的作用下就得到了图10(c)所示的采样器的输出信号 。图10 连续时间信号的采样过程梳状 函数 可以作为理想采样开关，来描述理想的采样过程。从数学上讲，采样信号 可以看成 和 的乘积，即: 从物理意义上，上述描述的理想采样过程可理解为脉冲调制过程：输入信号 作为调制信号，单位脉冲序列 作为载波，经过理想采样开关，得到理想采样序列 ，如图11所示。图11 理想采样过程3.2 采样定理信号的采样确定了连续时间信号 的采样表达式 ，那么，采样间隔 必须符合什么样的条件时， 才能保留有原连续时间信号 的所有信息。香农（Shannon）采样定理解决了此类问题。设连续时间信号 ，其傅里叶换为 ， 频谱中的最高频率成分为 ，对连续时间信号 采样，采样频率为 ，采样后的离散时间信号为 ，如果满足条件 ，可以从离散时间信号 中恢复原连续时间信号 ，否则，会发生频率混叠，从离散信号 中无法恢复原连续时间信号 。我们可以用下面的图解说明这个过程。假设 为连续函数， 为无穷脉冲序列，其采样间隔为 ，波形如图12(a), (b) 所示。将 与 相乘即可以得到离散的 信号，其波形如图12(e) 所示，此过程称为采样，也称为时域离散， 叫做采样函数。可以利用卷积定理来讨论采样在频域中的变化过程。 和 的傅里叶变换分别表示于图12(c) 和(d) 中。依据卷积定理，图12(f)中的图像应是图12(c)和(d) 的频域函数 的卷积。由于 是无限个 函数，所以根据函数卷积性质，只需将 移到每个 函数位置上即可，可见离散后的信号的频谱是一个周期函数，只需观察其中的一个周期即可，它与连续函数 的傅里叶变换相同。图12 波形采样过程 （普通采样频率） 图13 波形采样过程（较低采样频率）如在上例中，加大采样间隔 (即采样频率变小)，其结果如图13所示。可见由于增大了 ，频域采样函数 变得更密，因为频域脉冲函数的间隔减小，它们与频谱 的卷积就产生了波形的相互重叠，如图13(f)所示，函数的傅里叶变换由于采样引起的这种畸变称为混叠效应。混叠效应产生的原因在于采样间隔 太大， 也就是采样频率过低。如何才能不产生这种现象呢？在图14中可以看出，当 脉冲函数的频率间隔小于 ，即 ，卷积便出现重叠现象，这里 是连续函数 的最高频率成分。图14 不产生混叠现象时的频谱如果令 ， 称为采样频率， 不发生混叠效应的公式为： (15)即： 如果认为信号中的最高频率 等于 的话， 则有： (16)这就是香农采样定理。实际工程中，取： 对未知的信号最高频率成分 ，可由信号的低通滤波器的截止频率产生。 3.3 采样点数与频率分辨率信号分析处理过程中，首先要确定的是采样频率 ，采样频率应符合采样定理。其次要确定采样点数 和频率分辨率 。采样点数受到快速傅里叶变换( FFT) 变换的限制，一般情况下，只能取1024或2048 这类 的点数。采样频率 、采样点数 和频率分辨率 之间的关系如下： (17)式中， 为采样总时间长度（单位为s）； 为采样间隔。频率分辨率决定了频谱的分析精度，要想提高频率分辨率（降低 ），应尽量减少采样频率 ，或增加采样点数 。例，已知某信号的最高频率 ，希望达到的频率分辨率 ，试选择采样频率 、采样长度 及采样 。采样频率 ，取： 频率分辨率 ，采样时间长度 为： 采样点数 ，于是： 此时， 应取 应取1024，采样频率为 ： 由于采样长度未变，此时，频率分辨率仍然不变，且 符合要求。4、离散傅里叶变换（DFT） 离散傅里叶变换 (Discrete Fourier Transform, DFT)是一种用计算机计算傅里叶变换的方法。4.1 DFT的理论公式及计算过程连续时间信号 经采样后得到长度为 的时间列 ，此时频谱 仍是连续的，但是，实际运算中只能对有限项进行计算，因此必须对连续的频率轴离散化，以便与时域采样信号相对应。根据DFT的定义，一个连续函数的DFT求解过程需要三步：时域离散、时域截断和频域离散，时域离散的过程已经在3.1中讨论过了，这里仅简单介绍，重点讨论时域截断和频域离散后两步过程。(1) 第一步：时域离散如图15(a)~(c)所示，对 [图15(a)]进行波形采样，即将 乘以采样函数 [图15(b)]，采样间隔为 。图15(c)表示采样结果和它的傅里叶变换。这是对原始波形的第一次修改，一般将这一步修改叫做时域离散。由前面的3.2可知，采样后多少会产生一些混叠误差，可遵循采样定理来减少这一误差。(2) 第二步：时域截断采样后的函数 ，仍有无穷多个样本值，不适合计算机计算，所以必须将采样后的函数 进行截断，使之为有限个样本点。图15(d)表示了截断函数(矩形窗函数)和它的傅里叶变换。 为截断函数的区间宽度。图15(e)表示截断后有限宽度的离散时间函数，其傅里叶变换是带混叠效应的频率函数与截断窗函数的傅里叶变换作卷积。它的影响已在图15(e)中表示出来，使其傅里叶变换结果出现了皱褶(为了清楚起见，图中有意夸大了)，称为频率泄漏。要想减少这种频率泄漏，可将截断窗函数的长度尽可能选得长些，因为矩形函数的傅里叶变换是 ，若增加截断函数的长度 ，该函数就越近于 函数。这时，在频域的函数卷积中由于截断所引起的皱褶(或误差)也就越小。一般总希望截断函数的长度尽可能选得长些，但是，由于计算机内存的限制，不可能无限长，因此，频率泄漏现象总是存在的，这在后面我们还要讨论。这是对原信号傅里叶变换的第二次修改，称为对离散波形的时域截断，时域截断会产生频域泄漏误差。 图15 DFT过程的图解说明(3) 第三步：频域离散(或称时域延拓)将原始波形修改到图15(e)还不行，因为频域函数还是连续的，仍不能用计算机进行频域的计算。所以还必须用频域采样函数 把频域函数离散化。根据式(17)，频率离散的间隔 应为 。我们知道，时域采样引起频域的周期化。那么，在频域上采样的结果也会引起时域函数的周期化。当用图15(f)所示频域采样函数 与图15(e)所示的频域函数相乘时，相当于图15(f)所示时域函数 与 相卷积，由于 也是梳状 函数，相当于把 平移至各个 函数处，如图15(g)的 所示，这引起了时域函数的周期化，称为时域延拓。为了保证时域延拓时周期函数能够首尾相接，频域采样间隔 必须是窗长度 的倒数。频域采样也会引起一种误差，称为栅栏效应。我们可以把频域采样函数出现的位置想象成光栅的一根线，用其与原始信号相比，只有光栅线位置上的信号能被看到，其他都被栅栏挡住了，所以称为栅栏效应。与时域信号不同，频谱多呈函数的特性，所以很小的栅栏效应就会引起很大的幅值误差。解决栅栏效应的办法就是减少 的值，即减少采样频率，或加长时域截断窗 的长度，但这些也都是有限制的。因此，栅栏效应也总是多少存在的。图15(g)为最终的变换结果，可见，离散里叶变换相当于将原时间函数和频域函数二者都截断并修改成周期函数。其中中间的 个时间采样值和 个频域采样值为真值其余均为周期化产生的。或者说，为了进行傅里叶变换，DFT方法用长度为 矩形函数截取任意一段非周期函数，将其左右时域延拓后，形成一个新的周期为 的周期函数，然后对其进行傅里叶变换，根据傅里叶级数的特点，此时信号的基频 即为频谱信号的频率分辨率。值得一提的是，当采用计算机计算傅里叶变换时，得到的频谱并不在图15(g)所示的位置。而在图16所示的位置。当频谱的离散时间序列为 时，其中 为正频率部分，根据频谱的延拓性， 部分其实为负频率部分，或者说是右边第一个延拓谱的负频率部分。不管怎么说，由于实数的频谱是正负频率对称的，因此，只有前一半的计算点( 点)有意义。图16 DFT运算频谱的实际位置 5、快速傅里叶变换（FFT） 5.1 快速傅里叶变换(FFT)的方法提出背景FFT其是一种高效快速计算DFT的算法。对于 点的时间序列 的直接计算DFT，需要进行 次复数乘法和 复数加法的运算。已知一次复数乘法等于四次实数乘法，一次复数加法等于两次实数加法。因此，对大的 来说(数据处理中一般取 =1024)，这是一个相当大的运算量。所以，虽然早就有了DFT理论及计算方法，但因计算工作量大、计算时间长，限制了实际应用，迫使人们想办法提高对DFT的计算速度。1965年，美国学者库利(Cooly)和图基(Tueky)提出了快速算法，即FFT(Fast Fourier Transform，FFT)算法。基于FFT总的计算量为 次，相对直接DFT减少次数为 ，可见计算次数减少一半以上。5.2 FFT计算案例对一个时间函数 进行分析，其中信号频率 采样频率 ，采样点数 点，计算结果如图17所示 图17 FFT运算结果可见，在图17(b) 点频谱图上除了10Hz 处的谱峰外，还有 90 Hz 谱峰，这是10Hz 的负频率成分，即被测10Hz信号的幅值轴对称成分，因此 FFT 的运算结果只取前 点即可，即图 7(a) 所示的 点频谱图是正确的。## 导入包from matplotlib import pyplot as pltfrom matplotlib import rcParamsimport numpy as npimport pandas as pdconfig = { "font.family": 'serif', # 衬线字体 "font.size": 14, # 相当于小四大小 "font.serif": ['SimSun'], # 宋体 "mathtext.fontset": 'stix', # matplotlib渲染数学字体时使用的字体，和Times New Roman差别不大 'axes.unicode_minus': False # 处理负号，即-号}rcParams.update(config) ##========绘制时域信号图========##def plt_time_domain(arr, fs=1600, ylabel='Amp(mg)', title='原始数据时域图', img_save_path=None, x_vline=None, y_hline=None): """ :fun: 绘制时域图模板 :param arr: 输入一维数组数据 :param fs: 采样频率 :param ylabel: y轴标签 :param title: 图标题 :return: None """ import matplotlib.pyplot as plt plt.rcParams['font.sans-serif'] = ['SimHei'] # 显示中文 plt.rcParams['axes.unicode_minus'] = False # 显示负号 font = {'family': 'Times New Roman', 'size': '20', 'color': '0.5', 'weight': 'bold'} plt.figure(figsize=(12,4)) length = len(arr) t = np.linspace(0, length/fs, length) plt.plot(t, arr, c='g') plt.xlabel('t(s)') plt.ylabel(ylabel) plt.title(title) if x_vline: plt.vlines(x=x_vline, ymin=np.min(arr), ymax=np.max(arr), linestyle='--', colors='r') if y_hline: plt.hlines(y=0.2, xmin=np.min(t), xmax=np.max(t), linestyle=':', colors='y') #===保存图片====# if img_save_path: plt.savefig(img_save_path, dpi=500, bbox_inches = 'tight') plt.show() ##========绘制频域信号图========##def plt_fft_img(arr, fs, ylabel='Amp(mg)', title='频域图', img_save_path=None, vline=None, hline=None, xlim=None): """ :fun: 绘制频域图模板 :param arr: 输入一维时域数组数据 :param fs: 采样频率 :param ylabel: y轴标签 :param title: 图标题 :return: None """ # 计算频域幅值 length = len(arr) t = np.linspace(0, length/fs, length) fft_result = np.fft.fft(arr) fft_freq= np.fft.fftfreq(len(arr), d=t[1]-t[0]) # FFT频率 fft_amp= 2*np.abs(fft_result)/len(t) # FFT幅值 # 绘制频域图 plt.figure(figsize=(12,4)) ##=====单边频谱图====## plt.subplot(1,2,1) plt.plot(fft_freq[0: int(len(t)/2)], fft_amp[0: int(len(t)/2)], label='Frequency Spectrum', color='b') plt.xlabel('频率 (Hz)') plt.ylabel('幅值') plt.legend() ##=====双边频谱图====## plt.subplot(1,2,2) plt.plot(fft_freq[0: int(len(t))], fft_amp[0: int(len(t))], label='Frequency Spectrum', color='b') plt.xlabel('频率 (Hz)') plt.ylabel('幅值') plt.legend() if vline: plt.vlines(x=vline, ymin=np.min(fft_amp), ymax=np.max(fft_amp), linestyle='--', colors='r') if hline: plt.hlines(y=hline, xmin=np.min(fft_freq), xmax=np.max(fft_freq), linestyle=':', colors='y') #===保存图片====# if img_save_path: plt.savefig(img_save_path, dpi=500, bbox_inches = 'tight') if xlim: # 图片横坐标是否设置xlim plt.xlim(0, xlim) plt.tight_layout() plt.show() return fft_freq, fft_amp fs = 100 # 采样频率f = 10 # 模拟正弦信号频率time = 10.24 # 采样时长t = np.linspace(0, time, int(time*fs))data = 10*np.sin(2*np.pi*f*t) plt_time_domain(data, fs)fft_freq, fft_amp = plt_fft_img(data, int(time*fs)) 注意，幅值为6点多，与10相差了一些，这是因为时域截断时没有截到整周期，导致出现了频率泄露，因此幅值降低了。可以试试把采样时长设置为10s，这个时候幅值接近10，即频率泄露很轻了。 6、总结 频域是什么，频域包含哪些信息？在实际应用场景中，很多信号都为周期信号，如旋转机械振动、交流电流信号。这个周期其实也等同于频率，把这个周期信号用频率、幅值和相位进行描述，就是到频域了。 时域是怎么转换到频域的，时域转换到频域的原理是什么？傅里叶级数告诉我们，任何一个周期函数都可以用无限个三角函数去拟合。时域信号通过分解成傅里叶级数，即一个时域信号是由哪些个三角函数组合，把它给找出来，这个过程就是傅里叶变换。连续信号怎么转换到离散信号，有什么要注意的吗？因为计算机只能离散信号，所以需要将连续信号离散化。连续信号经过等周期的定期采样，得到离散信号。信号离散化需要满足采样定理，即采样频率 需要大于2倍信号中最大频率 。 FT、DFT、FFT之间是什么关系？FT，傅里叶变换，是一种将信号从时域转换到频域的理论，研究的主要是连续信号。DFT，离散傅里叶变换，由于工程应用需要，计算机只能处理离散信号，因此DFT主要是处理离散信号转换至频域。FFT，快速傅里叶变换，由于DFT直接计算比较复杂，计算工作量大，因此提出一种能提升计算效率的算法，即FFT。频率混叠、频率泄露分别是什么？当采样频率 小于2倍信号中最大频率 ，实际频谱中的高频部分会出现叠加，即频率混叠。在DFT时，需要用一个窗口截取信号，该窗口较小，或这窗口没有截取到周期部分时，导致在时域延拓时，信号首尾没有对接上，在频谱上出现了其它额外的频率成分，即出现了频率泄露。 直流分量是什么，会有什么影响？当信号中有直流部分，如 中1就是直流部分。当信号中含有直流成分时，在频谱中的0频率附近就会有较高幅值的谱峰出现，这会严重影响其他计算结果的显示比例。 编辑：李正平审核、校对：陈凯歌，赵栓栓，曹希铭，赵学功，白亮，陈少华如需转载，请后台联系作者说明：部分图片来源网络，若有侵权，烦请联系处理来源：故障诊断与python学习
