不用傅里叶变换，提取某一频率的幅值和相位

不用傅里叶变换，提取某一频率的幅值和相位

摘要： 本文从实际工程问题入手，探寻解决办法，为引入信号正交分解，和广义傅里叶级数做铺垫。

转子做周期性旋转时，不平衡质量所产生的周期性惯性离心力会引起转子产生简谐运动。在工程现场测得某一不平衡转子产生的振动加速度信号(记为$sig$)如下图所示：

从图中我们可以观测出不平衡分量的频率大约为7.4 Hz, 现在抛出问题，如何提取该信号不平衡成分的幅值和相位？
大家首先想到的肯定是：傅里叶变换

傅里叶变换肯定是可以的，但今天讨论的是，假如没有傅里叶变换，该如何解决？本文会给出求解过程并引出信号处理中的重要理论——信号的正交分解和广义傅里叶级数。求解过程如下：

假设由该不平衡量产生的振动加速度波形(记为$fit$)可以表示为：
$$fit=A\cdot cos(2\pi ft+\phi )$$
则该问题转换为在$f=7.4$ Hz时，寻找一对$\left(\begin{array}{c}A,\phi \end{array}\right)$使得，$fit$与上述图中的信号$sig$距离$\varphi$最小：

$$\varphi (A,\phi )=\sum _{j=1}^N(fit[j]-sig[j]{)}^2$$
其中$N$为振动信号$sig$和$fit$的点数，为方便计算，对fit进行展开，有
$$fit=A\cdot cos(\phi )\cdot cos(2\pi ft)-A\cdot sin(\phi )\cdot sin(2\pi ft)$$

令 $a=A\cdot cos(\phi )$ , $b=-A\cdot sin(\phi )$, $x=cos(2\pi ft)$, $y=sin(2\pi ft)$, 则有

$$fit=a\cdot x+b\cdot y$$
因为 $f$ 是确定的，采样时刻序列 $t$ 也是确定的，所以 $x$ 和 $y$ 都是常量序列。则此时 $fit$ 与 $sig$ 的距离 $ \phi $ 可以表示为以 $a$，$b$ 为变量的表达式 ：

$$\varphi (a,b)=\sum _{j=1}^N(fit[j]-sig[j]{)}^2$$

将 $fit=a\cdot x+b\cdot y$ 代入上式，有：
$$\varphi (a,b)=\sum _{j=1}^N(a\cdot x[j]+b\cdot y[j]-sig[j]{)}^2$$
将上式完全展开，有：

$$\varphi (a,b)=\sum _{j=1}^N(a^2\cdot x^2[j]+b^2\cdot y^2[j]+si{g}^2[j]+2abx[j]y[j]-2a\cdot x[j]sig[j]-2b\cdot y[j]sig[j])$$

根据连续函数极值定理，在 $\varphi$ 最小时，有：
$$\frac{\partial \varphi }{\partial a}=\frac{\partial \varphi }{\partial b}=0$$

整理，得到方程组：

$$\{\begin{array}{l}{\sum }_{j=1}^N{x}^2[j]\cdot a+{\sum }_{j=1}^{N}x[j]y[j]\cdot b={\sum }_{j=1}^{N}sig[j]x[j]\\ {\sum }_{j=1}^{N}y[j]x[j]\cdot a+{\sum }_{j=1}^N{y}^2[j]\cdot b={\sum }_{j=1}^{N}sig[j]y[j]\end{array}$$

令：

$$X=\left[\begin{array}{cc}x[1]& y[1]\\ x[2]& y[2]\\ ...& ...\\ x[N]& y[N]\end{array}\right],Y=\left[\begin{array}{c}sig[1]\\ sig[2]\\ ...\\ sig[N]\end{array}\right],\alpha =[a,b]$$

则上述方程组可表示为：

$$X^{T}X\alpha =X^{T}Y$$

进而求得：
$$\alpha =(X^{T}X{)}^{-1}{X}^{T}Y$$

求得 $\alpha$ 后，便得到了 $a,b$ 的值，进而通过下面的公式得到$A,\phi$ 的值：
$$A=\sqrt{a^2+b^2},\phi =-arctan\frac{b}{a}$$
至此，问题得到解决，上述方法成功提取到了该信号不平衡成分的幅值和相位。效果见下图所示，核心代码在文末给出。

这时，有人会问了，你是怎么观测到不平衡分量的频率大约为7.4 Hz的？我看着怎么有点像7.3Hz呢？其实相差不大，下图是本文分析过程分别在f为7.4, 7.3, 7.2Hz时的结果对比图。

那么如何理解这种结果呢？

其实，当f等于7.4Hz时，所提取的信号成分是原始信号在"7.4Hz方向"上的投影。7.4Hz方向和7.3Hz方向以及7.2Hz方向的"夹角"很小，所以投影的结果差距也不是很大，因此在精度要求不是特别苛刻时，可以使用肉眼观测的频率(7.4, 7.3, 7.2都可以)进行提取。

通过上面的论述可以知道，信号可以在任意多个方向上进行投影，在很多时候，我们需要的可能是信号在多个方向上的投影，那么有没有一种通用的大家都认可的投影方式呢？

答案是肯定的，参考二维平面中向量的表示方法或许可以得到一些启示：二维平面中的向量可以用一对坐标来表示，例如$（3,4）$ 代表该向量在x轴方向上投影的长度为3，在y轴方向上投影的长度为4. 那么该向量可以分解为x轴方向和y轴方向的两个向量$（3,0）$和$（0,4）$，在中学学习受力分析时我们也是这么干的。
那么为什么沿着x轴方向和y轴方向呢？因为这两个方向是正交的。而正交分解得到的分量具有很多美好的性质，在下一篇文章中，会深入浅出地介绍信号的正交分解和广义傅里叶级数。
微信公众号“振动信号研究所”会在第一时间发布新作，关注后回复“unbalance”获取本文所用波形数据。

import numpy as np
import matplotlib.pyplot as plt

sig = np.loadtxt(r".\unbalanced.txt")
Fs = 8192
t = np.arange(len(sig))/Fs
f = np.array([7.4 for _ in t])  # 欲提取的频率成分
x = np.cos(f*t*2*np.pi)
y = np.sin(f*t*2*np.pi)

X = np.mat(np.asarray([x, y]).T)
Y = np.mat(sig).T
xtx = X.T * X

alpha = xtx.I * (X.T * Y)

a = alpha.tolist()[0][0]
b = alpha.tolist()[1][0]
A = np.sqrt(a**2 + b**2)
phi = - np.arctan(b/a) if - np.arctan(b/a) < 0 else -np.pi - np.arctan(b/a)
fit = A*np.cos(2*np.pi*f*t + phi)
plt.figure()
plt.plot(t, sig, label="sig")
plt.plot(t, fit, label="fit")
plt.title(f"fit: {np.round(A, 3)}*np.cos(2*np.pi*{f[0]}*t+({np.round(phi, 3)}))")
plt.legend(loc="upper right")
plt.show()