36、MATLAB 中信号分析与频谱处理实践

MATLAB 中信号分析与频谱处理实践

1. 频率分析的初步策略

在进行信号频率分析时，若存在疑问，建议先在较宽的频率范围和较低的频率分辨率下进行分析。这有助于我们对信号的整体频率分布有一个初步的认识。基于此分析结果，我们可以进一步缩小感兴趣的频率范围的上下限。之后，增加频谱中的线条数量，从而评估更高频率分辨率对频谱的影响。

2. 测试信号与 MATLAB 实现

测试信号存储在相关材料中，主要以 AUDIO 文件形式存在，可通过常见的 PC 声卡播放。这些测试信号可用于测试 FFT 分析仪和测量仪器。几乎所有测试信号的采样频率为 44.1 kHz，只有一个测试信号的采样频率为 96 kHz，该信号由两个正弦波组成，频率分别为 f1 = 100 Hz 和 f2 = 45 kHz，振幅值为 0.354，用于测试混叠滤波器。

在 MATLAB 中，基于傅里叶变换将信号从时域转换到频域的语句如下：

% 此处为原文档中的语句 (14:78)，但文档未明确给出具体内容

结果 s 包含一个具有 N 个复数的向量。不过，实际执行的 DFT 没有乘以 1/N（归一化）。对于线数 N = 2^n（如 1024、2048 等），使用了计算优化算法。MATLAB 使用 FFTW [30] 程序库进行傅里叶变换。

为了获得可用的变换结果，需要进行归一化：

% 此处为原文档中的语句 (14:79)，但文档未明确给出具体内容

此时会得到一个双边复频谱，通过以下语句可以显示双边振幅谱：

% 此处为原文档中的语句 (14:80)，但文档未明确给出具体内容

MATLAB 函数不需要采样频率 fs 的信息，因此不会提供用于缩放横坐标的频率向量 f。通过以下语句可以呈现具有正确振幅值和频率轴的单边振幅谱：

% 此处为原文档中的语句 (14:81)，但文档未明确给出具体内容

帕塞瓦尔定理（时域功率等于频谱功率）仍然有效，但单边频谱中只有 N/2 条线，这导致对原公式 (14:54) 进行修正：

P = (1/N) * sum(x.^2) == 2 * sum(abs(X(1:N/2)).^2); % 对应原文档中的语句 (14:82)

由于单边频谱缺乏镜像，其幅值需乘以 2。

为了验证处理过程的正确性，对由九个具有倍频程中心频率的正弦波组成的测试信号进行 1/n 倍频程滤波，然后确定峰值。

3. 傅里叶变换函数的选择

如果想遵循某些建议，MATLAB 函数 fft() 较难使用。在实际应用中，更通用的傅里叶变换函数如下：

% 此处为原文档中的语句 (14:83)，但文档未明确给出具体内容

如果时间数据向量 x 比傅里叶变换所需的值的数量（参数 N）长，则将时间数据向量划分为多个长度为 N 的数据块，并分别进行傅里叶变换。在进行傅里叶变换之前，数据块需乘以窗函数 w，至少需要矩形窗。也可以为参数 w 指定整数值 N，其在分割中的效果与矩形窗函数相同。此时，s 变为包含单边复频谱的矩阵，p 包含功率谱密度矩阵，t 变为可用于缩放时间轴的向量，这种三维表示称为频谱图。

4. 频谱图的创建步骤

以测试信号 fs_44_1 kHz_Octaves.wav 为例，创建频谱图的步骤如下：
| 语句 | 解释 |
| ---- | ---- |
| df = 0.5 | 确定频率分辨率 Δf |
| N = ceil(fs/df) | 二分频谱的线数，使用 ceil() 函数取整 |
| T = N/fs | 块长度 |
| overlap = ceil(2/3 N) | 重叠部分，必须为整数值 |
| w = hann(N) | 窗函数，可在 MATLAB 帮助中查找其他窗函数 |
| FM = sum(w)/N | 根据公式 (14.68) 计算窗函数的平均值 |
| PM = sum(w.^2)/N | 根据公式 (14.64) 计算平均功率值 |
| B_eff = PM/(T * FM^2) | 根据公式 (14.65) 计算有效带宽 |
| [s, f, t, p] = spectrogram(x, w, overlap, N, fs) | 计算频谱图。若不需要 PSDs，可用波浪号 (~) 作为占位符 |
| sNorm = s/N | 单边复频谱的归一化 |
| MAG = 2 abs(sNorm)/FM | 计算振幅谱 |
| surf(t, f, MAG) | 显示振幅频谱图 |
| colormap(spring) | 定义颜色梯度，可在 MATLAB 帮助中查找其他梯度 |
| c = colorbar | 显示颜色条用于振幅值的颜色缩放 |
| set(gca, ‘FontSize’, 16) | 设置字体大小 |
| title(‘Spectrogram of test signal: octaves’) | 频谱图标题 |
| xlabel(‘time in s’) | X 轴标签 |
| ylabel(‘Frequency in Hz’) | Y 轴标签 |
| set(gca, ‘yscale’, ‘log’) | 将 Y 轴转换为对数刻度 |
| zlabel(‘Amplitude’) | Z 轴标签 |
| view(120, 20) | 旋转频谱图以便更好地可视化 |
| axis([0 30 0 20000 0 0.15]) | 定义轴缩放的值范围 |
| caxis([0 0.15]) | 调整颜色缩放以匹配轴缩放 |

5. 实际应用示例 - 柴油发动机振动测量

为了确定两缸柴油发动机的运行平稳性，使用加速度计进行振动测量，测量信号的采样频率 fs = 44100 Hz，选择的频率分辨率 Δf = 2 Hz，得到块长度 T = 0.5 s。其余傅里叶变换参数按照上述表格中的步骤选择。

通过频谱图可以观察到轻微的速度波动和测量过程中的速度变化。线性平均是快速了解信号中包含哪些频率的方法，从频谱图中可以读取到两个频率：
- 60 Hz：一阶发动机频率，由点火顺序和一阶自由质量力产生。
- 120 Hz：二阶发动机频率，由二阶自由质量力产生。

对于频率范围高达 1000 Hz 的情况，以不同视角显示的三维频谱图可以更清晰地展示频率峰值的产生原因。例如，从“上方”观察频谱图可以发现，频率峰值仅由最后 10 秒的激励引起，并且在测量过程中电机负载增加，电机速度在第 20 到 25 秒之间的变化比前后更明显。

6. 频率扫描信号的频谱分析

测试信号 fs_44_1 kHz_sweep.wav 是一个振幅为 1 的余弦信号，在 30 秒内频率从 f1 = 25 Hz 变化到 f2 = 2500 Hz。频谱图以三种不同的视图显示，包括从“上方”观察的整体频谱图，以及从右侧和频率轴方向观察的水平线条图。

在频率表示中，振幅线呈现梳状，这是由于频率分辨率不适合该扫描信号。该扫描信号的频率变化斜率为：

f_ST = (2500 - 25) / 30; % 对应原文档中的语句 (14:87)

当 DFT 的 Δf = 0.5 Hz 时，块长度 T = 2 s，在块长度内信号的频率变化为：

Delta_f_Signal = T * f_ST; % 对应原文档中的语句 (14:88)

这对应于频谱中的 330 条线。虽然 66.7% 的重叠可以稍微提高振幅值，但总体而言，参数设置不正确。为了获得真实的振幅值，Δf 至少应采用 f_ST = 82.5 Hz。

当 Δf = 82.5 Hz 时，DFT 的结果更接近信号的真实振幅值，但在 f 接近 0 和 t 接近 0 的范围内，振幅值过高，这是由于信号形状（余弦）、使用的窗函数（汉宁窗）、重叠操作（66.7%）以及由此产生的振幅校正的综合影响。解决方法是舍弃信号的起始部分。

7. 速度上升/下降过程的频谱分析

对于之前示例中的电机，在速度上升/下降过程中进行测量，持续时间为 180 秒，在前 150 秒内，速度从 n = 700/min 上升到 n = 4000/min。根据公式计算二阶电机频率的变化率 f_ST = 0.74 Hz/s，因此可以使用 Δf = 2 Hz 计算该信号的 DFT。

频谱图显示了速度上升和下降过程中的频率和时间切片。时间切片的索引通过向量 t 中的时间间隔确定，计算公式为：

% 假设 t 为时间向量，t1 为第一个频谱的标称时间，Delta_t 为时间间隔
Index = (t - t1) / Delta_t + 1; % 对应原文档中的语句 (14:91)

频谱平均可以提高结果的显著性。对于频谱图中的 1078 个频谱，时间间隔 ΔtSpek = 0.1666 s，根据频率变化率 f_ST 和选择的 DFT 参数 Δf = 2 Hz 以及 66.7% 的重叠，可以计算出在信号频率变化 2 Hz 时经过的频谱数量，从而进行无信息损失的平均。平均数量的计算公式为：

n = (Delta_f / f_ST - T) / Delta_tSpek + 1; % 对应原文档中的语句 (14:92)

取向下取整后的整数作为可用于平均的频谱数量。对于该频谱图，可能的平均数量 n = 14.2，平均时间 TA = 2.7 s，BeffTA 乘积为 4.05，相对标准偏差 ε 为 24.8%。在 MATLAB 代码中，平均操作通过循环实现：

% 此处为原文档中的语句 (14:94)，但文档未明确给出具体内容

除了平均频谱外，还需要形成新的时间向量 t2。如果最后一次平均使用的频谱少于 14 个，则按实际情况处理。

另一种平均方式是从第 7 个频谱开始，平滑地进行平均直到倒数第 7 个频谱，这种方式会在频谱图的开头和结尾损失 14 个频谱。MATLAB 代码如下：

% 此处为原文档中的语句 (14:95)，但文档未明确给出具体内容

还可以将测量时间划分为离散的时间间隔 tStep，并为每个时间点分配 n 个平均值作为频谱。通常将频谱的标称时间设置在块长度 T 的中间，在该时间点左右各取 n 个频谱进行平均。MATLAB 代码中相关设置如下：

% 此处为原文档中的语句 (14:96)，但文档未明确给出具体内容

根据这些设置确定逐步计算频谱图的参数：

% 此处为原文档中的语句 (14:97)，但文档未明确给出具体内容

通过循环在信号上按间隔 nStep 调用 spectrogram 指令计算频谱并进行平均，同时计算时间向量 t。

综上所述，通过合理选择频率分析策略、正确设置 MATLAB 中的参数以及采用适当的频谱平均方法，我们可以更准确地分析信号的频率特性，为实际应用提供有力的支持。

graph LR
    A[开始] --> B[初步频率分析]
    B --> C[缩小频率范围]
    C --> D[增加频谱线数]
    D --> E[选择测试信号]
    E --> F[MATLAB 时域到频域转换]
    F --> G[归一化处理]
    G --> H[显示频谱图]
    H --> I[分析信号频率特性]
    I --> J[判断是否需要平均]
    J -- 是 --> K[计算平均数量]
    K --> L[进行频谱平均]
    J -- 否 --> M[结束分析]
    L --> M

以上是关于信号分析和频谱处理的详细介绍，希望对大家有所帮助。在实际应用中，需要根据具体情况灵活调整参数和方法，以获得最佳的分析结果。

MATLAB 中信号分析与频谱处理实践

8. 频谱平均的不同策略对比

在信号分析中，频谱平均是提高结果显著性的重要手段。前文介绍了几种频谱平均的方法，下面对它们进行详细对比。

• 基于频率变化的平均法 ：根据信号的频率变化率 (f_{ST}) 和选择的频率分辨率 (\Delta f) 计算可用于平均的频谱数量。例如在速度上升/下降的电机测量中，通过公式 (n = \frac{\Delta f}{f_{ST}} - \frac{T}{\Delta t_{Spek}} + 1) 计算平均数量。这种方法适用于信号频率有明显变化的情况，能够在保证信息不损失的前提下提高结果的准确性。但需要准确测量信号的频率变化率，并且计算过程相对复杂。
• 从第 7 个频谱开始的平滑平均法 ：从第 7 个频谱开始，平滑地进行平均直到倒数第 7 个频谱。这种方法会在频谱图的开头和结尾损失一定数量的频谱，但能使平均过程更加平滑。适用于对频谱图开头和结尾部分数据要求不高，更注重整体平滑性的情况。
• 离散时间间隔平均法 ：将测量时间划分为离散的时间间隔 (t_{Step})，为每个时间点分配 (n) 个平均值作为频谱。通常将频谱的标称时间设置在块长度 (T) 的中间，在该时间点左右各取 (n) 个频谱进行平均。这种方法适用于需要对特定时间点的频谱进行详细分析的情况，能够更灵活地处理不同时间点的频谱信息。

平均方法	适用情况	优点	缺点
基于频率变化的平均法	信号频率有明显变化	信息不损失，结果准确	计算复杂，需准确测量频率变化率
从第 7 个频谱开始的平滑平均法	注重整体平滑性	平均过程平滑	开头和结尾损失部分频谱
离散时间间隔平均法	需要对特定时间点频谱详细分析	灵活处理不同时间点频谱信息	计算量较大

9. 窗函数的选择与影响

在进行傅里叶变换之前，数据块需要乘以窗函数 (w)。窗函数的选择会对频谱分析结果产生重要影响。

• 矩形窗 ：是最简单的窗函数，至少在进行傅里叶变换时需要使用。为参数 (w) 指定整数值 (N) 的效果与矩形窗函数相同。矩形窗的特点是主瓣窄，频率分辨率高，但旁瓣高，会产生频谱泄漏。
• 汉宁窗（Hanning） ：在示例中使用的窗函数，其表达式为 (w = hann(N))。汉宁窗的主瓣比矩形窗宽，频率分辨率相对较低，但旁瓣低，能够有效抑制频谱泄漏。适用于大多数信号分析场景，尤其是对频谱泄漏较为敏感的情况。

不同窗函数的选择应根据具体的信号特点和分析要求来决定。如果需要高频率分辨率，可以选择矩形窗；如果需要抑制频谱泄漏，汉宁窗是更好的选择。

10. 信号分析的误差来源与解决方法

在信号分析过程中，可能会出现各种误差，影响分析结果的准确性。下面分析常见的误差来源及解决方法。

• 频率分辨率不合适 ：如在频率扫描信号分析中，当频率分辨率 (\Delta f) 不适合信号的频率变化时，会导致振幅线呈现梳状，无法准确反映信号的真实振幅。解决方法是根据信号的频率变化率 (f_{ST}) 合理选择频率分辨率，使 (\Delta f) 至少不小于 (f_{ST})。
• 频谱泄漏 ：由于窗函数的影响，信号的频谱会泄漏到相邻的频率范围，导致频谱分析结果不准确。选择合适的窗函数可以有效抑制频谱泄漏，如汉宁窗。
• 平均参数设置不当 ：在频谱平均过程中，如果平均数量 (n) 或平均时间 (T_A) 设置不当，会影响结果的显著性。需要根据信号的特点和分析要求，合理计算平均数量和平均时间。

11. 总结与实践建议

通过对信号分析和频谱处理的详细介绍，我们可以总结出以下实践建议：

• 频率分析策略 ：在进行信号频率分析时，先在宽频率范围和低频率分辨率下进行初步分析，然后逐步缩小频率范围，增加频谱线数，以获得更准确的频谱信息。
• MATLAB 函数选择 ：根据具体需求选择合适的傅里叶变换函数，对于复杂的分析任务，更通用的函数可能更合适。
• 参数设置 ：合理设置频率分辨率、窗函数、重叠率等参数，根据信号的特点和分析要求进行调整。
• 频谱平均 ：根据信号的频率变化情况和分析目的，选择合适的频谱平均方法，提高结果的显著性。

在实际应用中，需要不断尝试和调整参数，结合具体情况进行分析，以获得最佳的分析结果。

graph LR
    A[选择窗函数] --> B{窗函数类型}
    B -- 矩形窗 --> C[高频率分辨率，高旁瓣]
    B -- 汉宁窗 --> D[低旁瓣，抑制频谱泄漏]
    E[分析误差来源] --> F{误差类型}
    F -- 频率分辨率不合适 --> G[调整频率分辨率]
    F -- 频谱泄漏 --> H[选择合适窗函数]
    F -- 平均参数设置不当 --> I[重新计算平均参数]
    J[总结实践建议] --> K[合理频率分析策略]
    J --> L[选择合适 MATLAB 函数]
    J --> M[合理设置参数]
    J --> N[选择合适频谱平均方法]

总之，信号分析和频谱处理是一个复杂而重要的领域，通过合理运用 MATLAB 工具和相关方法，可以深入了解信号的频率特性，为实际应用提供有力的支持。希望本文的介绍能够帮助读者更好地进行信号分析和频谱处理。