37、信号分析：频域处理、调制分析与传递函数

信号分析：频域处理、调制分析与传递函数

1. 信号频域分析基础

在信号处理中，对于无时间变化的信号，常进行平均处理。频谱平均的频谱数量考量，可通过特定公式（Eq. 14.68）实现。除线性平均外，还有其他频谱平均或评估方法，如下表所示：

功能	函数	描述
线性均值	mean(s, dim)	按指定计算方向 dim 计算矩阵 s 的线性均值，公式为 (x = \frac{1}{n}\sum_{\xi = 1}^{n}x(\xi))
二次均值	rms(s, dim)	按指定计算方向 dim 计算矩阵 s 的均方根值，公式为 (X_{RMS} = \sqrt{\frac{1}{n}\sum_{\xi = 1}^{n}x(\xi)^2})
总和	sum(s, dim)	按指定计算方向 dim 计算矩阵 s 的总和
标准差	std(s, [], dim)	按指定计算方向 dim 计算矩阵 s 的标准差
最大值	max(s, [], dim)	按指定计算方向 dim 确定矩阵 s 中的最大值
众数	mode(s, dim)	按指定计算方向 dim 确定矩阵 s 中最常出现的值

MATLAB 函数计算的频谱矩阵是按列排列的，即第一列存储第一个频谱，第二列存储第二个频谱，依此类推。因此，频谱的平均或其他评估必须逐行进行，这可通过在相应指令中添加参数 dim = 2 来实现。

1.1 频谱量的 MATLAB 实现

频谱中的纵坐标值会根据信号类型、使用场景和约定而有所不同。MATLAB 代码（Table 14.13）的实现并非简单的纵坐标重新缩放。通过特定的 MATLAB 语句，可将信号 x 从时域转换到频域，得到归一化的单边复谱图 sNorm。在特殊情况下，若时间信号的值数量与所需线数 N/2 相同，则谱图可简化为频谱，后续无需进行平均等处理。

1.1.1 幅度谱

线性平均幅度谱可通过特定指令计算，结果为平均峰值 (x(f))。对于均方根值，可通过相应公式计算，并以 dBV 为单位显示为电平。为实现所需的幅度校正，需将归一化频谱 sNorm 乘以 1/FM（FM 依据 Eq. 14.63）。

1.1.2 功率谱

平均功率谱可根据 MATLAB 指令计算。为使结果更易理解和操作，可将指令分步执行：
1. 从归一化复频谱矩阵 s 形成幅度谱。
2. 计算归一化复频谱矩阵的均方根值。
3. 将均方根值除以 (\sqrt{2})。
4. 对结果进行平方得到功率谱，再通过线性平均简化为频谱。为正确表示功率，幅度谱需乘以校正值 1/FM。

1.1.3 功率谱密度（PSD）

确定 PSD 有三种方法：
1. 通过 MATLAB 函数 spectrogram 并进行后续线性平均。
2. 使用 MATLAB 函数 pwelch，该函数已包含平均功能。
3. 按 Eq. 14.68 分步执行指令。

这三种计算方法的结果相同。

1.1.4 倒谱

倒谱可通过 MATLAB 函数计算。

2. 调制分析

幅度调制是将低频振荡叠加到高频振荡上，导致振幅随时间变化。这种调制振荡常见于机器中，如齿轮箱，可能由轴偏心等因素触发。

2.1 调制分析步骤

进行调制分析需以下步骤：
1. 使用 1/n 倍频程滤波将信号划分为不同的频带。
2. 为每个频带形成包络。
3. 对包络进行离散傅里叶变换（DFT）。

最终结果以彩色缩放的二维图像呈现，横坐标为载波频率，纵坐标为调制频率，颜色信息表示调制深度。

2.2 示例：齿轮箱调制分析

以一个齿轮级为例，包含齿数为 40 的齿轮 Z1 和齿数为 25 的齿轮 Z2。当 Z1 的转速为 3000/min（旋转频率 (f_{D,Z1} = 50 Hz)）时，齿轮啮合频率 (f_Z = 2000 Hz)；Z2 的转速为 4800/min（旋转频率 (f_{D,Z2} = 80 Hz)）。通过幅度调制，齿轮可对齿啮合频率进行调制。从调制分析结果可知，载波频率约为 2000 Hz 时，受到约 50 Hz 的强烈调制，可推断出齿轮啮合频率 (f_Z = 2000 Hz) 受齿轮 Z1 的旋转频率（(f_{D,Z1} = 50 Hz)）调制。

2.3 MATLAB 代码实现

为避免在后续傅里叶变换中分析滤波器瞬态，需将滤波后的信号和时间向量在开头和结尾各缩短 5 s。对于 1/3 倍频程滤波，需定义滤波器参数。在循环中，依次进行信号滤波、包络计算和傅里叶变换，形成每个频带的调制幅度谱。最终，变量 MOD 包含一个 6401 × 24 的矩阵，代表调制频率与载波频率带的关系。通过相应指令，可将调制分析结果进行图形显示，并对横坐标的载波频率带标注正确的频率值。

2.4 示例：电动执行器调制分析

以车辆座椅调节电动执行器的空气声记录为例进行调制分析，结果显示多个频带受到调制频率 (f_{MOD} = 101 Hz) 的调制。为确定坐标轴的范围，可先熟悉显示结果，再逐步调整缩放比例，直至得到最终显示效果。

3. 传递函数

传递函数用于关联振荡系统的输入和输出变量，可完整描述线性时不变系统（LTI 系统），并据此推断系统本身及其参数。

3.1 传递函数的表示

在时域中，若输入 (x(t)) 和输出 (y(t)) 以时间函数形式给出，则有 (y(t) = \int_{-\infty}^{\infty}x(t)h(t - \tau)d\tau)，其中 (h(t)) 称为脉冲响应函数，符号表示为 (y(t) = x(t) * h(t))。但在测量实践中，通常更倾向于使用频域表示，即 (H(j\omega) = \frac{Y(j\omega)}{X(j\omega)})，其中 (H(j\omega)) 称为频率响应函数（FRF）。

3.2 传递函数的实验确定

以摩托车前叉的垂直振动为例，输入 (x(t)) 为轮胎接触点的激励，输出 (y(t)) 为振动响应。同步采集输入和输出信号，进行傅里叶变换得到复频谱 (X(j\omega)) 和 (Y(j\omega))，进而形成传递函数 (H(j\omega))。由于叠加噪声会扭曲传递函数，因此需对传递函数进行修正。通过将分子和分母乘以复共轭 (X^ )，得到传递函数 (H_1)：
[H_1(j\omega) = \frac{Y(j\omega) \cdot X^ (j\omega)}{X(j\omega) \cdot X^*(j\omega)} = \frac{S_{XY}(j\omega)}{S_{XX}(\omega)} = \frac{G_{XY}(j\omega)}{G_{XX}(\omega)}]

同理，将分子和分母乘以复共轭 (Y^ )，可得到传递函数 (H_2)：
[H_2(j\omega) = \frac{Y(j\omega) \cdot Y^ (j\omega)}{X(j\omega) \cdot Y^*(j\omega)} = \frac{S_{YY}(\omega)}{S_{YX}(j\omega)} = \frac{G_{YY}(\omega)}{G_{YX}(j\omega)}]

在信号分析仪的用户界面或软件中，需选择使用 (H_1) 或 (H_2)。

3.3 噪声对传递函数的影响

当振荡系统的输入或输出叠加干扰信号时，(H_1) 和 (H_2) 的差异会显现出来，具体情况如下：
- 输出有噪声 ：系统受 (e(t)) 激励，输入测量信号 (x(t)) 无干扰，输出信号 (a(t)) 叠加干扰 (m(t)) 或 (n(t))，且噪声与激励无关。此时，传递函数 (H_2) 可提供无系统误差的传递函数 (H) 的估计。
- 输入有噪声 ：激励 (e(t)) 叠加干扰 (m(t))，输入测量信号 (x(t)) 有噪声，输出信号 (a(t)) 无干扰。此时，传递函数 (H_1) 可提供仅含统计误差的传递函数 (H) 的估计，而 (H_2) 包含系统误差 (G_{nn}/G_{aa})。
- 输入和输出都有噪声 ：(H_1) 和 (H_2) 都存在系统误差，(H_1) 返回的值偏小，(H_2) 返回的值偏大，传递函数 (H) 介于 (H_1) 和 (H_2) 之间。通常，(H_1) 对最小值的近似效果更好，(H_2) 对最大值的近似效果更好。

3.4 相干函数

为评估输出频谱 (Y) 相对于输入频谱 (X) 的线性相关性，可使用相干函数 (\gamma^2)：
[\gamma^2(\omega) = \frac{|S_{YX}(j\omega)|^2}{S_{XX}(\omega) \cdot S_{YY}(\omega)}]

相干函数 (\gamma^2) 为实数，取值范围在 0 到 1 之间。其值表示在该频率下，信号 (y(t)) 源自信号 (x(t)) 的程度，可用于判断两个信号之间的因果关系。相干函数与实验传递函数 (H_1) 和 (H_2) 之间存在关系 (\gamma^2 = \frac{H_1(j\omega)}{H_2(j\omega)})。

当 (\gamma^2 = 1) 时，信号完全线性相关，系统为线性时不变系统；当 (\gamma^2 = 0) 时，信号完全线性无关，因果关系不存在。实际测量中，常得到 (0 < \gamma^2 < 1) 的值，原因可能包括输入和/或输出叠加干扰、信号之间的非线性关系、窗函数导致的泄漏以及信号之间的传输时间差异等。通过优化测量链，可使相干函数在感兴趣的频率范围内最大化。相干函数还可解释为相应频率下的信噪比 (SNR = \frac{\gamma^2}{1 - \gamma^2})。

3.5 传递函数的应用示例

3.5.1 动态质量的确定

以摆锤为例，用脉冲锤敲击悬挂在两根线上的紧凑质量，通过测量加速度和力，将测量信号转换到频域，形成传递函数 (H_1(f))。若将加速度作为输入（单位：(m/s^2)），力作为输出（单位：N），则动态质量 (m_{dyn}) 可通过 (m_{dyn} = H(f) = \frac{F}{\ddot{x}}) 计算。由于测量是在紧凑质量上进行的，预计幅度频率响应具有频率独立性。

为进行信号分析，需将输入和输出信号截断，以确保包含整个脉冲锤敲击过程。由于脉冲信号在窗口边界处值为零，可使用矩形窗口，无需进行频谱校正。测量结果显示，在 20 Hz 至约 5 kHz 的频率范围内，传递函数的幅度频率响应几乎恒定；在高于 5 kHz 的频率范围内，由于加速度和力频谱的弯曲，传递函数逐渐趋近于“零除以零”的形式，导致结果出现偏差。相干函数在整个频率范围内值为 1，表明线性相关性，但在高于 5 kHz 的频率范围内，传递函数的结果明显错误。这提示在实际应用中，激励必须在整个感兴趣的频率范围内具有足够的功率，以获得足够的响应信号。

3.5.2 飞轮共振频率的确定

通过传递函数可确定飞轮的声音决定共振频率及其阻尼。具体过程与摆锤测量类似，使用 MATLAB 函数 modalfit 可精确确定共振频率。从测量结果可知，飞轮的共振频率约为 1222 Hz、3215 Hz 和 3235 Hz。

4. 传递函数相关操作流程总结

4.1 传递函数实验测定流程

下面通过 mermaid 流程图展示传递函数实验测定的详细流程：

graph TD
    A[同步采集输入 x(t) 和输出 y(t) 信号] --> B[对 x(t) 和 y(t) 进行傅里叶变换]
    B --> C[得到复频谱 X(jω) 和 Y(jω)]
    C --> D{是否有噪声干扰}
    D -- 无噪声 --> E[直接计算 H(jω)=Y(jω)/X(jω)]
    D -- 有噪声 --> F{噪声位置}
    F -- 输出有噪声 --> G[计算 H2(jω)=GYY(ω)/GYX(jω)]
    F -- 输入有噪声 --> H[计算 H1(jω)=GXY(jω)/GXX(ω)]
    F -- 输入输出都有噪声 --> I[H 介于 H1 和 H2 之间，综合考量]

4.2 调制分析操作流程

调制分析的操作步骤可以总结为以下表格：
| 步骤 | 操作内容 |
| — | — |
| 1 | 使用 1/n 倍频程滤波将信号划分为频带 |
| 2 | 为每个频带形成包络 |
| 3 | 对包络进行离散傅里叶变换（DFT） |
| 4 | 以彩色缩放二维图像呈现结果，横坐标为载波频率，纵坐标为调制频率，颜色表示调制深度 |

5. 不同频谱计算方法对比

5.1 功率谱密度（PSD）计算方法对比

确定功率谱密度（PSD）有三种方法，下面对这三种方法进行详细对比：
| 方法 | 具体操作 | 优点 | 缺点 |
| — | — | — | — |
| 方法一 | 通过 MATLAB 函数 spectrogram 并进行后续线性平均 | 步骤清晰，易于理解 | 计算过程相对复杂，需要额外进行线性平均 |
| 方法二 | 使用 MATLAB 函数 pwelch，该函数已包含平均功能 | 计算简便，函数自带平均功能 | 对函数内部实现细节依赖较大 |
| 方法三 | 按 Eq. 14.68 分步执行指令 | 灵活性高，可以根据具体需求调整步骤 | 操作步骤较多，容易出错 |

5.2 不同频谱类型计算要点

不同类型的频谱计算有各自的要点，总结如下：
- 幅度谱 ：线性平均幅度谱计算结果为平均峰值 (x(f))，均方根值计算后以 dBV 为单位显示电平，需将归一化频谱乘以 1/FM 进行幅度校正。
- 功率谱 ：分步计算，先形成幅度谱，再计算均方根值，除以 (\sqrt{2}) 后平方得到功率谱，最后线性平均简化为频谱，幅度谱需乘以校正值 1/FM 以正确表示功率。
- 功率谱密度（PSD） ：三种计算方法结果相同，可根据实际需求选择合适的方法。
- 倒谱 ：通过 MATLAB 函数计算。

6. 实际应用中的注意事项

6.1 传递函数测量中的噪声处理

在传递函数测量中，噪声是一个不可忽视的问题。当输入或输出叠加噪声时，会影响传递函数的准确性。为了减少噪声的影响，应尽量消除干扰源，确保测量信号 (x(t)) 和 (y(t)) 尽可能无中间干扰源。如果无法避免噪声，应优化测量设置，将干扰源分配到输出或输入端。同时，结合相干函数进行分析，当相干函数值较低时，传递函数的结果可能不可靠，此时应谨慎解读。

6.2 频谱计算中的参数选择

在频谱计算中，参数的选择至关重要。例如，在使用 MATLAB 函数进行频谱计算时，矩阵的计算方向参数 dim 需要根据频谱矩阵的排列方式正确设置。对于幅度谱、功率谱等的计算，校正因子 1/FM 的选择要依据具体的公式（如 Eq. 14.63）。在调制分析中，1/n 倍频程滤波的参数设置会影响频带的划分，进而影响后续的分析结果。

6.3 实验激励的功率要求

从传递函数的实际应用案例中可以看出，激励必须在整个感兴趣的频率范围内具有足够的功率，以获得足够的响应信号。特别是在存在阻尼和非线性因素（如 backlash、friction、非线性弹簧和阻尼特性等）的情况下，更需要保证激励功率，否则传递函数的结果可能会出现偏差。

7. 总结

本文详细介绍了信号分析中的频域处理、调制分析和传递函数相关知识。在频域处理方面，阐述了频谱平均的方法、不同频谱类型（幅度谱、功率谱、功率谱密度、倒谱）的计算方法及 MATLAB 实现。调制分析部分介绍了幅度调制的原理、调制分析的步骤及实际应用案例。传递函数部分讲解了其定义、实验确定方法、噪声对其的影响、相干函数的作用以及实际应用示例。

通过对这些内容的学习，我们可以更好地理解信号的特性，解决实际工程中的信号分析问题。在实际应用中，要注意噪声处理、参数选择和激励功率等问题，结合相干函数等指标综合分析传递函数的结果，以获得准确可靠的信号分析结论。同时，合理运用 MATLAB 等工具进行频谱计算和分析，能够提高工作效率和分析精度。