最后更新: 2026-03-10 07:28 查看: 111 次纠错评论(0) 课件开VIP

信号的提取

> 在[上文](https://kb.kmath.cn/kbase/detail.aspx?id=3174)说过，一个函数可以投影到正交三角函数系上，就可以获得该分量上的函数强度，如果把函数看成信号，这句话也可以理解为信号在这个分量上的强度，简称滤波。

## 不用傅里叶变换采样
转子做周期性旋转时，不平衡质量所产生的周期性惯性离心力会引起转子产生简谐运动。在工程现场测得某一不平衡转子产生的振动加速度信号(记为$ sig$)如下图所示：
  ![图片](/uploads/2025-07/17a528.jpg)
  
 
图中横坐标代表时间，从图中肉眼可以观察到$1s$内振动大约7.4次，即频率大约为 $f=7.4 Hz$ ,（周期为频率的倒数所以周期$T=0.13s$），现在抛出问题，**如何提取该信号周期性信号的的幅值和相位?**
 
 大家首先想到的肯定是：傅里叶变换

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

假设由该不平衡量产生的振动加速度波形（记为 $f i t$ ）可以表示为：

$$
f i t=A \cdot \cos (2 \pi f t+\varphi)
$$

则该问题转换为在 $f=7.4 Hz$ 时，寻找一对 $(A, \varphi)$ 使得，fit与上述图中的信号 sig距离 $\phi$ 最小，即

$$
\phi(A, \varphi)=\sum_{j=1}^N(f i t[j]-s i g[j])^2
$$

其中 $N$ 为振动信号 $\operatorname{sig}$ 和 $f i t$ 的点数，为方便计算，对fit进行展开，有
$$
f i t=A \cdot \cos (\varphi) \cdot \cos (2 \pi f t)-A \cdot \sin (\varphi) \cdot \sin (2 \pi f t)
$$

令 $a=A \cdot \cos (\varphi), b=-A \cdot \sin (\varphi), x=\cos (2 \pi f t), y=\sin (2 \pi f t)$ ，则有

$$
f i t=a \cdot x+b \cdot y
$$

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

$$
\phi(a, b)=\sum_{j=1}^N(f i t[j]-\operatorname{sig}[j])^2
$$

将 $f i t=a \cdot x+b \cdot y$ 代入上式，有：

$$
\phi(a, b)=\sum_{j=1}^N(a \cdot x[j]+b \cdot y[j]-s i g[j])^2
$$

将上式完全展开，有：

$$
\begin{gathered}
\phi(a, b)=\sum_{j=1}^N\left(a^2 \cdot x^2[j]+b^2 \cdot y^2[j]+\operatorname{sig}^2[j]+2 a b x[j] y[j]-2 a \cdot x[j] \operatorname{sig}[j]\right. -2 b \cdot y[j] \operatorname{sig}[j])
\end{gathered}
$$

根据连续函数极值定理，在 $\phi$ 最小时，有：

$$
\frac{\partial \dot{\phi}}{\partial a}=\frac{\partial \phi}{\partial b}=0
$$

整理，得到方程组：

$$
\left\{\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 \operatorname{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 \operatorname{sig}[j] y[j]
\end{array}\right.
$$

令：

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

则上述方程组可表示为：

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

进而求得：

$$
\alpha=\left(X^T X\right)^{-1} X^T Y
$$

求得 $\alpha$ 后，便得到了 $a, b$ 的值，进而通过下面的公式得到 $A, \varphi$ 的值：

$$
A=\sqrt{a^2+b^2}, \varphi=-\arctan \frac{b}{a}
$$

至此，问题得到解决，上述方法成功提取到了该信号不平衡成分的幅值和相位。效果见下图所示，

![图片](/uploads/2025-07/dfb864.jpg)
 
 
 这时，有人会问了，你是怎么观测到不平衡分量的频率大约为7.4 Hz的？我看着怎么有点像7.3Hz呢？其实相差不大，下图是本文分析过程分别在f为7.4, 7.3, 7.2Hz时的结果对比图。
 
 ![图片](/uploads/2025-07/4ee05a.jpg)
 
  
### 那么如何理解这种结果呢？
其实，当f等于 7.4 Hz 时，所提取的信号成分是原始信号在＂ 7.4 Hz 方向＂上的投影。 7.4 Hz 方向和 7.3 Hz 方向以及 7.2 Hz 方向的＂夹角＂很小，所以投影的结果差距也不是很大，因此在精度要求不是特别苛刻时，可以使用肉眼观测的频率（7．4，7．3，7．2都可以）进行提取。

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

答案是肯定的，参考二维平面中向量的表示方法或许可以得到一些启示：二维平面中的向量可以用一对坐标来表示，例如 $(3,4)$ 代表该向量在 $x$ 轴方向上投影的长度为 3 ，在 $y$ 轴方向上投影的长度为 4 ．那么该向量可以分解为 $x$ 轴方向和 $y$ 轴方向的两个向量 $(3,0)$ 和 $(0,4)$ ，在中学学习受力分析时我们也是这么干的。那么为什么沿着x轴方向和y轴方向呢？因为这两个方向是正交的。而正交分解得到的分量具有很多美好的性质，下面再介绍一下信号的正交分解和广义傅里叶级数。

## 信号的正交分解与广义傅里叶级数

早期电视的天线接收信号如下图

![图片](/uploads/2026-02/08919a.jpg){width=300px}

当信号通过天线时，会被天线截取，形成电磁信号。
 ![天线](/uploads/2026-02/25af9a.gif){width=400px}
 
如果把电磁信号进行抽样，变成离散的，就是信号提取，如下图
 ![电磁信号](/uploads/2026-02/ceddea.gif){width=400px}
 
 
假设有 $n$ 个函数 $\varphi_1(t), \varphi_2(t), \ldots \varphi_n(t)$ 在区间 $\left(t_1, t_2\right)$ 构成一个正交函数集 $\left\{\varphi_i(t)\right\}$ ，将任一函数 $f(t)$ 用这 $n$ 个正交函数的线性组合来近似，可以表示为：

$$
f(t) \approx C_1 \varphi_1(t)+C_2 \varphi_2(t)+C_n \varphi_n(t)=\sum_{j=1}^n C_j \varphi_j(t)
$$

估计量的误差定义为：

$$
e(t)=f(t)-\sum_{j=1}^n C_j \varphi_j(t)
$$

问题是如何选择各系数 $C_j$ 使 $f(t)$ 与近似函数之间的误差在区间 $\left(t_1, t_2\right)$ 内最小？通常，采用最小均方**误差准则**。估计量的均方误差定义为

$$
J[e(t)]=\frac{1}{t_2-t_1} \int_{t_1}^{t_2}\left|f(t)-\sum_{j=1}^n C_j \varphi_j(t)\right|^2 d t
$$

为使上式最小，由连续函数极值定理，有

$$
\frac{\partial J[e(t)]}{\partial C_i}=\frac{\partial}{\partial C_i}\left\{\frac{1}{t_2-t_1} \int_{t_1}^{t_2}\left|f(t)-\sum_{j=1}^n C_j \varphi_j(t)\right|^2 d t\right\}=0
$$

展开上式中的被积函数，并求导。结合正交函数集的性质，上式非 0 项为：

$$
\frac{\partial}{\partial C_i}\left\{\frac{1}{t_2-t_1} \int_{t_1}^{t_2}\left[-2 C_i \varphi_i^*(t) f(t)+C_i^2 \varphi_i(t) \varphi_i^*(t)\right] d t\right\}=0
$$

即

$$
-2 \int_{t_1}^{t_2} f(t) \varphi_i^*(t) d t+2 C_i \int_{t_1}^{t_2} \varphi_i(t) \varphi_i^*(t) d t=0
$$

所以，欲使均方误差最小，其第 $i$ 项函数 $\varphi_i(t)$ 的系数应该按照下式取值

$$
C_i=\frac{\int_{t_1}^{t_2} f(t) \varphi_i^*(t) d t}{\int_{t_1}^{t_2} \varphi_i(t) \varphi_i^*(t) d t}=\frac{1}{K_i} \int_{t_1}^{t_2} f(t) \varphi_i^*(t) d t
$$

在用正交函数去近似 $f(t)$ 时，所取得的项数越多，即 $n$ 越大，则均方误差越小。在 $n \rightarrow \infty$ 时（完备正交函数集），$f(t)$ 可以精确地用这 $n$ 个正交函数的线性组合来表示，如下式，此时均方误差为 0 ．

$$
f(t)=\lim _{n \rightarrow \infty} \sum_{i=1}^n C_i \varphi_i(t)
$$

这称为函数 $f(t)$ 的正交展开式，亦称为**广义傅里叶级数**。其中 $C_i=\frac{\int_{t_1}^{t_2} f(t) \varphi_i^*(t) d t}{\int_{t_1}^{t_2} \varphi_i(t) \varphi_i^*(t) d t}$ 称为**傅里叶系数**。

## 实践分析

我们用上面中的数据来实践本文所介绍的知识，上面文章欲从信号中提取的频率成分为 $f=7.4 Hz$ ，若要用本文所述方法来提取，则可以构造完备正交的三角函数集 $\{1, \cos (\Omega t), \cos (2 \Omega t), \ldots, \sin (\Omega t), \sin (2 \Omega t), \ldots\}$ ，将原信号（这里将信号看成连续函数）在区间 $(0, T)$ 上用所构造正交函数集的线性组合来表示，得到 7.4 Hz 对应的正、余弦函数前的系数便可以完成提取。这便要求

$$
\Omega=\frac{2 \pi f}{m}, \text { 其中: } m \in Z, f=7.4
$$

因为所用数据采样时长为 1 秒，为了利用到尽可能多的数据来求系数，则应在 $T=2 \pi / \Omega<1$ 的约束下，求 $T$ 的最大值。

求得当 $m=7$ 时， T 最大为 0.94595 ，则关系的系数为函数集中 $\cos (7 \Omega t), \sin (7 \Omega t)$ 对应的系数，分别记为 $a 、 b$ ，用本文方法求得：

$$
\begin{aligned}
a & =\frac{\int_0^{0.94595} \operatorname{sig}(t) \cos (7 \Omega t) d t}{\int_0^{0.94595} \cos ^2(7 \Omega t) d t} \approx-0.129 \\
b & =\frac{\int_0^{0.94595} \operatorname{sig}(t) \sin (7 \Omega t) d t}{\int_0^{0.94595} \sin ^2(7 \Omega t) d t} \approx 2.507
\end{aligned}
$$
将这两个同频成分合成余弦形式并展示，如下图：

![图片](/uploads/2025-07/b9eb07.jpg)

那么和上面文章相比，哪一个方法的提取效果更好呢？可以在评论区内留下你的看法。
  
  
  
 
 原文地址 微信公众号：[振动信号研究所](https://mp.weixin.qq.com/s/O2ft2ELNXshEI2wK3DZ-0g) 或 [知乎](https://zhuanlan.zhihu.com/p/550490748) 以及 [信号提取](https://zhuanlan.zhihu.com/p/553408889)
 
 
附录：代码

``` 
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()
```

科数题库(单机版) 会议室预约系统(book)

今日还可看 0 篇未注册用户每天查看4篇，注册用户每天8篇，开通VIP会员无限制查看。

免费注册

《高等数学》难点解析

高数教程 泰勒公式切线与法线切平面与法平面驻点·拐点·极值点·零点间断点渐进线瑕积分欧拉方程伯努利方程 Abel 收敛定理偏导数的几何意义偏导数的几何意义梯度数量场与向量场多元函数极值拉格朗日算子通量与散度环流量与旋度格林公式高斯公式斯托克斯公式三大公式比较傅里叶级数极坐标微元点法式方程变上限定积分 X型计算面积 Y型计算面积微分的意义渐近线间断点 y''+py'+qy=f(x)方程高斯黎曼傅里叶变换(复数) 拉普拉斯变换(复数)

《线性代数》难点解析

线代教程 近世代数对数学的整体思考线性的意义矩阵乘法(列视角) 矩阵乘法(行视角) 矩阵左乘矩阵右乘逆矩阵求解方程组阶梯形矩阵的求法方程组解的判定四阶行列式的计算线性变换的意义线性空间向量组的等价线性空间的几何意义基础解系的求法施密特正交化特征值与特征向量的意义矩阵相似的几何意义矩阵可对角化的理解秩的意义(向量版) 秩的意义(方程版) 二次型的意义

《概率论与数理统计》难点解析

概率教程 置信区间与上a分位数概率中的“取”与“放” 贝叶斯公式全概率公式泊松分布指数分布伽玛分布二维密度图的意义卷积的意义相关系数的意义 k阶矩是与矩母函数卡方分布的作用单正态区间估计理解假设检验理解切比雪夫不等式中心极限定理

上一篇：傅里叶变换例题

下一篇：傅里叶背景5：δ函数★★★★★

本文对您是否有用？有用 (0) 无用 (0)