【嵌入式 AI 实战第 5 期】工业振动检测（一）振动数据采集与信号分析

一、为什么工业振动检测是预测性维护的核心？

在现代制造业中，旋转机械设备（电机、泵、风机、齿轮箱、轴承）占工厂总设备数量的 60% 以上，其运行状态直接影响生产线的连续性和产品质量。据统计，80% 以上的旋转机械故障都会表现为振动异常，而传统的 "坏了再修" 和 "定期检修" 模式存在明显缺陷：

• 人工巡检效率低：依赖巡检人员经验，漏检率高达 30%，且无法 24 小时连续监测
• 故障发现滞后：往往在设备已经产生严重损坏时才被发现，维修成本增加 5-10 倍
• 过度维护浪费：定期更换完好部件，造成不必要的停机和备件浪费
• 安全隐患大：突发故障可能导致生产事故，甚至危及人员安全

预测性维护（PdM） 正是解决这些问题的关键技术，而振动检测则是预测性维护中应用最广泛、最成熟的技术手段。通过实时监测设备的振动信号，我们可以在故障萌芽阶段就发现异常，提前安排维修计划，实现 "按需维护"。

本期我们将从零开始，基于STM32F103C8T6+MPU6050搭建一套低成本、高性能的工业振动数据采集系统，深入讲解振动信号的预处理、滤波降噪和频谱分析方法。特别地，我们将提供完整的 C 语言实现，将所有信号处理算法直接运行在 STM32 嵌入式端，这才是工业现场的标准做法。最后我们将制作出可用于模型训练的工业振动故障数据集。

二、振动检测的核心理论基础

2.1 旋转机械振动产生机理

旋转机械的振动本质上是周期性的力作用在结构上产生的响应。当设备出现故障时，会产生特定频率的振动信号，这些信号就像设备的 "心电图"，包含了丰富的故障信息。

常见的故障类型及其振动特征：

2.2 MPU6050 惯性测量单元原理

MPU6050 是一款集成了三轴加速度计和三轴陀螺仪的 6 轴运动处理传感器，非常适合用于工业振动检测：

• 加速度计测量范围：±2g/±4g/±8g/±16g（工业振动检测推荐 ±8g）
• 输出数据速率：最高 1kHz（满足绝大多数旋转机械的采样需求）
• 内置 16 位 ADC，分辨率高
• 支持 I2C 通信，接口简单，易于与 STM32 连接

关键参数选择：

• 采样频率：根据奈奎斯特采样定理，采样频率应至少为被测信号最高频率的 2 倍。对于大多数工业电机（转速 3000rpm=50Hz），采样频率设置为 1kHz 即可覆盖前 20 次谐波
• 量程选择：工业设备正常振动加速度一般在 0.1-5g 范围内，选择 ±8g 量程可以兼顾灵敏度和动态范围

2.3 振动信号分析方法概述

振动信号分析主要分为时域分析和频域分析两大类：

时域分析：直接对时间序列信号进行分析，提取统计特征

• 基本统计量：均值、方差、标准差、峰值、峰峰值
• 无量纲指标：峭度、裕度、脉冲因子、波形因子
• 优点：计算简单，实时性好，对冲击性故障敏感

频域分析：通过傅里叶变换将时域信号转换为频域信号，分析不同频率成分的幅值

• 快速傅里叶变换（FFT）：最常用的频域分析方法
• 功率谱密度（PSD）：反映信号功率随频率的分布
• 优点：能够准确识别故障特征频率，定位故障部位

三、硬件系统设计与实现

3.1 系统总体架构

我们的振动数据采集系统支持两种工作模式：

传感器层(MPU6050) → 微控制器层(STM32F103) → 数据传输层(USB转串口)
↓
嵌入式DSP处理单元
↓
特征提取与报警逻辑

3.2 硬件电路设计

核心电路原理图

STM32F103C8T6 MPU6050
—————- ———-
PB6(I2C1_SCL) <–> SCL
PB7(I2C1_SDA) <–> SDA
3.3V <–> VCC
GND <–> GND
PA9(USART1_TX) <–> USB转串口RX
PA10(USART1_RX) <–> USB转串口TX
PA0 <–> 报警LED
PB0 <–> 蜂鸣器

工业现场硬件优化技巧

这是很多教程都会忽略的部分，但却是决定系统能否在工业现场稳定运行的关键：

3.3 硬件实物图与安装示例

这里展示实际搭建的硬件系统和正确的安装方式：

• 核心板：STM32F103C8T6 最小系统板
• 传感器：MPU6050 模块（带稳压电路）
• 通信：CH340 USB 转串口模块
• 报警：LED 指示灯 + 有源蜂鸣器
• 安装：使用热熔胶或螺丝固定在电机端盖上

四、STM32 基础软件实现

4.1 开发环境与工程配置

• 开发工具：STM32CubeMX 6.8.1 + Keil MDK 5.38
• 固件库：STM32CubeF1 1.8.4
• 配置步骤：
• 配置系统时钟为 72MHz
• 配置 I2C1（PB6、PB7），速率 100kHz
• 配置 USART1（PA9、PA10），波特率 115200
• 配置 TIM2 定时器，1ms 中断一次（用于精确采样）
• 配置 GPIO（PA0、PB0）用于报警输出
• 生成初始化代码

4.2 MPU6050 驱动实现

// mpu6050.h
#ifndef __MPU6050_H
#define __MPU6050_H

#include "stm32f1xx_hal.h"

#define MPU6050_ADDR 0xD0 // MPU6050 I2C地址

// 寄存器地址定义
#define PWR_MGMT_1 0x6B
#define SMPLRT_DIV 0x19
#define CONFIG 0x1A
#define GYRO_CONFIG 0x1B
#define ACCEL_CONFIG 0x1C
#define ACCEL_XOUT_H 0x3B
#define ACCEL_YOUT_H 0x3D
#define ACCEL_ZOUT_H 0x3F

// 加速度计量程配置
#define ACCEL_RANGE_2G 0x00
#define ACCEL_RANGE_4G 0x08
#define ACCEL_RANGE_8G 0x10
#define ACCEL_RANGE_16G 0x18

// 灵敏度系数
#define ACCEL_SENS_2G 16384.0f
#define ACCEL_SENS_4G 8192.0f
#define ACCEL_SENS_8G 4096.0f
#define ACCEL_SENS_16G 2048.0f

typedef struct {
int16_t accel_x_raw;
int16_t accel_y_raw;
int16_t accel_z_raw;
float accel_x; // 单位：g
float accel_y;
float accel_z;
} MPU6050_DataTypeDef;

uint8_t MPU6050_Init(I2C_HandleTypeDef *hi2c);
uint8_t MPU6050_ReadAccel(MPU6050_DataTypeDef *data);

#endif

// mpu6050.c
#include "mpu6050.h"

static I2C_HandleTypeDef *hi2c_mpu;
static float accel_sens;

uint8_t MPU6050_Init(I2C_HandleTypeDef *hi2c) {
uint8_t reg_val;
hi2c_mpu = hi2c;

// 唤醒MPU6050
reg_val = 0x00;
if (HAL_I2C_Mem_Write(hi2c_mpu, MPU6050_ADDR, PWR_MGMT_1, 1, &reg_val, 1, 1000) != HAL_OK) {
return 1;
}
HAL_Delay(100);

// 设置采样率分频器，采样率=1kHz/(1+分频器值)
reg_val = 0x00; // 采样率1kHz
HAL_I2C_Mem_Write(hi2c_mpu, MPU6050_ADDR, SMPLRT_DIV, 1, &reg_val, 1, 1000);

// 配置数字低通滤波器
reg_val = 0x03; // 截止频率44Hz
HAL_I2C_Mem_Write(hi2c_mpu, MPU6050_ADDR, CONFIG, 1, &reg_val, 1, 1000);

// 配置加速度计量程
reg_val = ACCEL_RANGE_8G;
HAL_I2C_Mem_Write(hi2c_mpu, MPU6050_ADDR, ACCEL_CONFIG, 1, &reg_val, 1, 1000);
accel_sens = ACCEL_SENS_8G;

return 0;
}

uint8_t MPU6050_ReadAccel(MPU6050_DataTypeDef *data) {
uint8_t buf[6];

if (HAL_I2C_Mem_Read(hi2c_mpu, MPU6050_ADDR, ACCEL_XOUT_H, 1, buf, 6, 1000) != HAL_OK) {
return 1;
}

data->accel_x_raw = (int16_t)((buf[0] << 8) | buf[1]);
data->accel_y_raw = (int16_t)((buf[2] << 8) | buf[3]);
data->accel_z_raw = (int16_t)((buf[4] << 8) | buf[5]);

// 转换为g单位
data->accel_x = data->accel_x_raw / accel_sens;
data->accel_y = data->accel_y_raw / accel_sens;
data->accel_z = data->accel_z_raw / accel_sens;

return 0;
}

五、上位机数据处理与可视化

5.1 Python 环境配置

我们使用 Python 进行数据的接收、处理和可视化，需要安装以下库：

pip install pyserial numpy matplotlib scipy pandas

5.2 串口数据接收与实时波形显示

import serial
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

# 串口配置
SERIAL_PORT = 'COM3' # 根据实际情况修改
BAUD_RATE = 115200
BUFFER_SIZE = 1000 # 显示的历史数据点数

# 初始化串口
ser = serial.Serial(SERIAL_PORT, BAUD_RATE, timeout=1)

# 初始化数据缓冲区
time_data = np.linspace(-BUFFER_SIZE/1000, 0, BUFFER_SIZE)
accel_x_data = np.zeros(BUFFER_SIZE)
accel_y_data = np.zeros(BUFFER_SIZE)
accel_z_data = np.zeros(BUFFER_SIZE)

# 创建图形
fig, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(12, 8), sharex=True)
line1, = ax1.plot(time_data, accel_x_data, 'r-', label='X轴加速度')
line2, = ax2.plot(time_data, accel_y_data, 'g-', label='Y轴加速度')
line3, = ax3.plot(time_data, accel_z_data, 'b-', label='Z轴加速度')

ax1.set_ylabel('加速度 (g)')
ax1.legend()
ax1.grid(True)
ax2.set_ylabel('加速度 (g)')
ax2.legend()
ax2.grid(True)
ax3.set_xlabel('时间 (s)')
ax3.set_ylabel('加速度 (g)')
ax3.legend()
ax3.grid(True)

plt.suptitle('实时振动加速度波形', fontsize=16)

def update(frame):
global accel_x_data, accel_y_data, accel_z_data

# 读取串口数据
while ser.in_waiting > 0:
line = ser.readline().decode('utf-8').strip()
try:
x, y, z = map(float, line.split(','))

# 更新数据缓冲区
accel_x_data = np.roll(accel_x_data, -1)
accel_y_data = np.roll(accel_y_data, -1)
accel_z_data = np.roll(accel_z_data, -1)

accel_x_data[-1] = x
accel_y_data[-1] = y
accel_z_data[-1] = z
except:
pass

# 更新图形
line1.set_ydata(accel_x_data)
line2.set_ydata(accel_y_data)
line3.set_ydata(accel_z_data)

# 自动调整Y轴范围
ax1.relim()
ax1.autoscale_view()
ax2.relim()
ax2.autoscale_view()
ax3.relim()
ax3.autoscale_view()

return line1, line2, line3

# 启动动画
ani = FuncAnimation(fig, update, interval=50, blit=True)
plt.tight_layout()
plt.show()

# 关闭串口
ser.close()

5.3 振动信号预处理算法

原始振动信号中包含大量的噪声和干扰，必须进行预处理才能用于后续的特征提取和分析。我们将实现三种常用的预处理算法：

1. 滑动窗口分割

将连续的时间序列信号分割成固定长度的窗口，每个窗口作为一个独立的样本进行处理。

def sliding_window(data, window_size, step_size):
"""
滑动窗口分割
:param data: 输入数据，形状为(n_samples, n_channels)
:param window_size: 窗口大小
:param step_size: 步长
:return: 分割后的数据，形状为(n_windows, window_size, n_channels)
"""
n_samples, n_channels = data.shape
n_windows = (n_samples – window_size) // step_size + 1

windows = np.zeros((n_windows, window_size, n_channels))

for i in range(n_windows):
start = i * step_size
end = start + window_size
windows[i] = data[start:end]

return windows

2. 去趋势滤波

去除信号中的直流分量和线性趋势，这些趋势通常是由传感器漂移或设备缓慢运动引起的。

from scipy.signal import detrend

def detrend_signal(data):
"""
去趋势滤波
:param data: 输入数据，形状为(n_windows, window_size, n_channels)
:return: 去趋势后的数据
"""
detrended_data = np.zeros_like(data)

for i in range(data.shape[0]):
for j in range(data.shape[2]):
detrended_data[i, :, j] = detrend(data[i, :, j])

return detrended_data

3. 均值滤波

简单有效的平滑滤波方法，能够去除高频噪声。

def mean_filter(data, kernel_size=5):
"""
均值滤波
:param data: 输入数据，形状为(n_windows, window_size, n_channels)
:param kernel_size: 卷积核大小
:return: 滤波后的数据
"""
filtered_data = np.zeros_like(data)
kernel = np.ones(kernel_size) / kernel_size

for i in range(data.shape[0]):
for j in range(data.shape[2]):
filtered_data[i, :, j] = np.convolve(data[i, :, j], kernel, mode='same')

return filtered_data

5.4 FFT 频谱分析

FFT 快速傅里叶变换是振动信号分析中最重要的工具，它能够将时域信号转换为频域信号，揭示信号的频率组成。

from scipy.fft import fft, fftfreq

def compute_fft(signal, sample_rate=1000):
"""
计算信号的FFT
:param signal: 输入时域信号，形状为(window_size,)
:param sample_rate: 采样频率，单位Hz
:return: 频率数组和幅值数组
"""
n = len(signal)
yf = fft(signal)
xf = fftfreq(n, 1 / sample_rate)[:n//2]
yf = 2.0 / n * np.abs(yf[:n//2])

return xf, yf

# 示例：对一个窗口的数据进行FFT分析
window_size = 1024
sample_rate = 1000

# 假设我们有一个窗口的X轴加速度数据
signal = detrended_data[0, :, 0]

# 计算FFT
xf, yf = compute_fft(signal, sample_rate)

# 绘制频谱图
plt.figure(figsize=(12, 6))
plt.plot(xf, yf)
plt.xlabel('频率 (Hz)')
plt.ylabel('幅值 (g)')
plt.title('振动信号频谱图')
plt.grid(True)
plt.xlim(0, 200) # 只显示0-200Hz的频率成分
plt.show()

注意：虽然上位机处理方便调试和数据分析，但在实际工业现场，我们通常会将信号处理算法直接在嵌入式设备端实现，这样可以获得更好的实时性和可靠性。具体实现方法见本文第八部分。

六、实验结果与分析

6.1 实验设置

我们使用一台普通的三相异步电机作为实验对象，设置了三种运行状态：

每种状态下采集 10 分钟的振动数据，采样频率 1kHz。

6.2 时域波形对比

从时域波形上可以明显看出三种状态的差异：

• 正常状态：波形平稳，幅值较小，一般在 ±0.2g 范围内
• 轻微故障状态：波形出现周期性波动，幅值增大到 ±0.5g 左右
• 严重故障状态：波形波动剧烈，幅值超过 ±1g，且有明显的冲击成分

6.3 频域频谱对比

频谱分析能够更清晰地揭示故障特征：

• 正常状态：频谱中主要是电机的转频（50Hz）及其低次谐波，幅值较小
• 轻微故障状态：转频（50Hz）的幅值明显增大，同时出现了 2 倍、3 倍转频
• 严重故障状态：转频及其谐波的幅值急剧增大，频谱中出现了丰富的高频成分

七、工业振动故障数据集制作

高质量的数据集是训练异常检测模型的基础。我们将按照以下步骤制作工业振动故障数据集：

7.1 数据集结构

vibration_dataset/
├── normal/
│ ├── sample_001.npy
│ ├── sample_002.npy
│ └── …
├── minor_fault/
│ ├── sample_001.npy
│ ├── sample_002.npy
│ └── …
└── severe_fault/
├── sample_001.npy
├── sample_002.npy
└── …

7.2 数据集制作脚本

import os
import numpy as np
import pandas as pd

# 配置参数
WINDOW_SIZE = 1024
STEP_SIZE = 512
SAMPLE_RATE = 1000

# 数据保存路径
DATASET_PATH = 'vibration_dataset'
os.makedirs(DATASET_PATH, exist_ok=True)
os.makedirs(os.path.join(DATASET_PATH, 'normal'), exist_ok=True)
os.makedirs(os.path.join(DATASET_PATH, 'minor_fault'), exist_ok=True)
os.makedirs(os.path.join(DATASET_PATH, 'severe_fault'), exist_ok=True)

def process_csv_file(csv_path, label):
"""
处理单个CSV文件，生成样本并保存
:param csv_path: CSV文件路径
:param label: 标签，'normal'、'minor_fault'或'severe_fault'
"""
# 读取CSV文件
df = pd.read_csv(csv_path, header=None, names=['x', 'y', 'z'])
data = df.values

# 滑动窗口分割
windows = sliding_window(data, WINDOW_SIZE, STEP_SIZE)

# 预处理
windows = detrend_signal(windows)
windows = mean_filter(windows)

# 保存每个窗口为单独的npy文件
save_dir = os.path.join(DATASET_PATH, label)
for i in range(windows.shape[0]):
save_path = os.path.join(save_dir, f'sample_{i:03d}.npy')
np.save(save_path, windows[i])

print(f"处理完成：{csv_path}，生成{windows.shape[0]}个样本")

# 示例：处理三个状态的数据
process_csv_file('normal_data.csv', 'normal')
process_csv_file('minor_fault_data.csv', 'minor_fault')
process_csv_file('severe_fault_data.csv', 'severe_fault')

# 生成数据集说明文件
with open(os.path.join(DATASET_PATH, 'README.md'), 'w') as f:
f.write('# 工业振动故障数据集\\n\\n')
f.write('## 数据集说明\\n')
f.write('- 采样频率：1000Hz\\n')
f.write(f'- 窗口大小：{WINDOW_SIZE}\\n')
f.write(f'- 步长：{STEP_SIZE}\\n')
f.write('- 传感器：MPU6050\\n')
f.write('- 量程：±8g\\n\\n')
f.write('## 标签说明\\n')
f.write('- normal：正常状态\\n')
f.write('- minor_fault：轻微故障状态\\n')
f.write('- severe_fault：严重故障状态\\n\\n')
f.write('## 数据格式\\n')
f.write('每个样本是一个形状为(1024, 3)的numpy数组，三列分别对应X、Y、Z轴加速度\\n')

print("数据集制作完成！")

7.3 数据集质量评估

制作完成后，我们需要对数据集进行质量评估：

八、嵌入式端信号处理算法 C 语言实现（工业现场标准方案）

在嵌入式设备端直接实现信号处理算法，相比上位机处理有三大核心优势：

• 实时性：无需传输原始数据，处理延迟从百毫秒级降至毫秒级
• 低带宽：只传输特征值而非原始数据，带宽占用降低 90% 以上
• 可靠性：不依赖网络和上位机，设备可独立运行和报警

下面我们提供完整的 C 语言实现，针对 STM32F103C8T6 优化，包含所有核心算法：滑动窗口、去趋势滤波、均值滤波、FFT 频谱分析和特征提取。同时提供CMSIS-DSP 库优化版本（工业项目首选）和纯 C 基础版本（无依赖）。

8.1 嵌入式端整体架构设计

我们采用采集 – 处理 – 输出流水线架构，使用环形缓冲区管理数据，实现边采集边处理，无数据丢失。

定时器中断采集(1ms) → 环形缓冲区 → 预处理(去趋势+均值滤波) → FFT频谱分析 → 特征提取 → 串口输出特征值

关键参数配置（针对 STM32F103C8T6 优化）：

• 采样频率：1000Hz
• FFT 点数：512 点（平衡精度和计算量，1024 点在 F1 上也能跑但稍慢）
• 滑动窗口步长：256 点（50% 重叠，提高数据利用率）
• 内存占用：约 12KB（完全在 F103 的 20KB SRAM 范围内）

8.2 基础数据结构与环形缓冲区

首先实现一个高效的环形缓冲区，用于存储连续采集的振动数据。

// vibration_dsp.h
#ifndef __VIBRATION_DSP_H
#define __VIBRATION_DSP_H

#include "stm32f1xx_hal.h"
#include <math.h>

#define SAMPLE_RATE 1000.0f // 采样频率(Hz)
#define FFT_POINTS 512 // FFT点数(必须是2的幂)
#define BUFFER_SIZE 1024 // 环形缓冲区大小(至少2倍FFT点数)

// 三轴加速度数据结构
typedef struct {
float x;
float y;
float z;
} AccelData;

// 环形缓冲区结构
typedef struct {
AccelData data[BUFFER_SIZE];
uint16_t head;
uint16_t tail;
uint16_t count;
} RingBuffer;

// 时域特征结构
typedef struct {
float mean; // 均值
float std; // 标准差
float peak; // 峰值
float rms; // 均方根
float kurtosis; // 峭度(冲击故障敏感指标)
float crest_factor;// 峰值因子
} TimeDomainFeatures;

// 频域特征结构
typedef struct {
float peak_freq; // 峰值频率
float peak_amp; // 峰值幅值
float total_power; // 总能量
float centroid; // 频谱质心
} FreqDomainFeatures;

// 环形缓冲区操作函数
void RingBuffer_Init(RingBuffer *rb);
uint8_t RingBuffer_Push(RingBuffer *rb, AccelData *data);
uint8_t RingBuffer_ReadWindow(RingBuffer *rb, AccelData *window, uint16_t window_size);

// 信号处理函数
void detrend(float *signal, uint16_t length);
void mean_filter(float *signal, uint16_t length, uint8_t kernel_size);
void compute_fft_magnitude(float *input, float *output, uint16_t n);
void arm_fft_magnitude(float *input, float *output, uint16_t n);

// 特征提取函数
void extract_time_domain_features(float *signal, uint16_t length, TimeDomainFeatures *features);
void extract_freq_domain_features(float *fft_mag, uint16_t n, FreqDomainFeatures *features);

#endif

// vibration_dsp.c
#include "vibration_dsp.h"

// 初始化环形缓冲区
void RingBuffer_Init(RingBuffer *rb) {
rb->head = 0;
rb->tail = 0;
rb->count = 0;
}

// 向缓冲区写入一个数据点
uint8_t RingBuffer_Push(RingBuffer *rb, AccelData *data) {
if (rb->count >= BUFFER_SIZE) {
return 1; // 缓冲区满
}

rb->data[rb->head] = *data;
rb->head = (rb->head + 1) % BUFFER_SIZE;
rb->count++;

return 0;
}

// 从缓冲区读取一个窗口的数据(先进先出)
uint8_t RingBuffer_ReadWindow(RingBuffer *rb, AccelData *window, uint16_t window_size) {
if (rb->count < window_size) {
return 1; // 数据不足
}

for (uint16_t i = 0; i < window_size; i++) {
window[i] = rb->data[(rb->tail + i) % BUFFER_SIZE];
}

// 移动尾指针(步长=window_size/2，实现50%重叠)
rb->tail = (rb->tail + window_size/2) % BUFFER_SIZE;
rb->count -= window_size/2;

return 0;
}

8.3 核心信号处理算法 C 实现

1. 去趋势滤波

去除传感器漂移和设备缓慢运动引起的线性趋势。

// 对单通道信号进行去趋势处理
void detrend(float *signal, uint16_t length) {
float sum_x = 0.0f, sum_y = 0.0f, sum_xy = 0.0f, sum_x2 = 0.0f;

// 计算线性拟合参数 y = a*x + b
for (uint16_t i = 0; i < length; i++) {
sum_x += i;
sum_y += signal[i];
sum_xy += i * signal[i];
sum_x2 += i * i;
}

float a = (length * sum_xy – sum_x * sum_y) / (length * sum_x2 – sum_x * sum_x);
float b = (sum_y – a * sum_x) / length;

// 减去线性趋势
for (uint16_t i = 0; i < length; i++) {
signal[i] -= (a * i + b);
}
}

2. 均值滤波

高效的滑动窗口均值滤波，去除高频随机噪声。

// 对单通道信号进行均值滤波
void mean_filter(float *signal, uint16_t length, uint8_t kernel_size) {
float sum = 0.0f;
uint8_t half_kernel = kernel_size / 2;

// 计算前kernel_size个点的和
for (uint8_t i = 0; i < kernel_size; i++) {
sum += signal[i];
}

// 滑动窗口计算
for (uint16_t i = half_kernel; i < length – half_kernel; i++) {
signal[i] = sum / kernel_size;
sum += signal[i + half_kernel + 1] – signal[i – half_kernel];
}
}

3. 纯 C 实现 FFT（无依赖版本）

这是一个基础的基 2 快速傅里叶变换实现，无需任何外部库，适合理解原理和资源极度受限的场景。

// 复数结构
typedef struct {
float real;
float imag;
} Complex;

// 复数乘法
static Complex complex_mul(Complex a, Complex b) {
Complex res;
res.real = a.real * b.real – a.imag * b.imag;
res.imag = a.real * b.imag + a.imag * b.real;
return res;
}

// 复数加法
static Complex complex_add(Complex a, Complex b) {
Complex res;
res.real = a.real + b.real;
res.imag = a.imag + b.imag;
return res;
}

// 复数减法
static Complex complex_sub(Complex a, Complex b) {
Complex res;
res.real = a.real – b.real;
res.imag = a.imag – b.imag;
return res;
}

// 位反转置换
static void bit_reverse(Complex *x, uint16_t n) {
uint16_t i, j, k;
for (i = 1, j = 0; i < n; i++) {
k = n >> 1;
while (j >= k) {
j -= k;
k >>= 1;
}
j += k;
if (i < j) {
Complex temp = x[i];
x[i] = x[j];
x[j] = temp;
}
}
}

// 基2FFT实现
void fft(Complex *x, uint16_t n) {
bit_reverse(x, n);

for (uint16_t s = 1; (1 << s) <= n; s++) {
uint16_t m = 1 << s;
float theta = -2 * M_PI / m;
Complex w_m = {cosf(theta), sinf(theta)};

for (uint16_t k = 0; k < n; k += m) {
Complex w = {1.0f, 0.0f};
for (uint16_t j = 0; j < m/2; j++) {
Complex t = complex_mul(w, x[k + j + m/2]);
Complex u = x[k + j];
x[k + j] = complex_add(u, t);
x[k + j + m/2] = complex_sub(u, t);
w = complex_mul(w, w_m);
}
}
}
}

// 计算FFT幅值谱
void compute_fft_magnitude(float *input, float *output, uint16_t n) {
Complex x[n];

// 初始化复数数组
for (uint16_t i = 0; i < n; i++) {
x[i].real = input[i];
x[i].imag = 0.0f;
}

// 执行FFT
fft(x, n);

// 计算幅值(只取前半部分)
for (uint16_t i = 0; i < n/2; i++) {
output[i] = 2.0f / n * sqrtf(x[i].real * x[i].real + x[i].imag * x[i].imag);
}
}

4. CMSIS-DSP 库优化版 FFT（工业项目首选）

强烈推荐在实际项目中使用 ARM 官方的 CMSIS-DSP 库，它针对 ARM 内核做了深度汇编优化，512 点 FFT 在 STM32F103 上仅需约12ms，比纯 C 实现快 5-10 倍。

第一步：添加 CMSIS-DSP 库到工程

第二步：优化版 FFT 实现

#include "arm_math.h"

// CMSIS-DSP优化版FFT
void arm_fft_magnitude(float *input, float *output, uint16_t n) {
arm_cfft_instance_f32 S;

// 初始化FFT实例
switch(n) {
case 128: arm_cfft_init_f32(&S, 128); break;
case 256: arm_cfft_init_f32(&S, 256); break;
case 512: arm_cfft_init_f32(&S, 512); break;
case 1024: arm_cfft_init_f32(&S, 1024); break;
default: return;
}

// CMSIS-DSP要求输入为复数数组(实部+虚部交替)
float fft_input[n*2];
for (uint16_t i = 0; i < n; i++) {
fft_input[2*i] = input[i];
fft_input[2*i+1] = 0.0f;
}

// 执行FFT
arm_cfft_f32(&S, fft_input, 0, 1);

// 计算幅值
arm_cmplx_mag_f32(fft_input, output, n);

// 归一化(只取前半部分)
for (uint16_t i = 0; i < n/2; i++) {
output[i] = 2.0f / n * output[i];
}
}

8.4 时域与频域特征提取

提取工业振动检测中最常用的特征指标，这些特征将直接用于下期的异常检测模型。

// 提取时域特征
void extract_time_domain_features(float *signal, uint16_t length, TimeDomainFeatures *features) {
float sum = 0.0f, sum_sq = 0.0f, sum_4th = 0.0f;
float max_val = -1e9f, min_val = 1e9f;

// 计算基本统计量
for (uint16_t i = 0; i < length; i++) {
sum += signal[i];
sum_sq += signal[i] * signal[i];
sum_4th += powf(signal[i], 4);

if (signal[i] > max_val) max_val = signal[i];
if (signal[i] < min_val) min_val = signal[i];
}

features->mean = sum / length;
features->rms = sqrtf(sum_sq / length);
features->peak = fmaxf(fabsf(max_val), fabsf(min_val));
features->std = sqrtf((sum_sq – sum*sum/length) / (length-1));

// 计算无量纲指标
features->crest_factor = features->peak / features->rms;
features->kurtosis = (sum_4th / length) / powf(features->rms, 4) – 3.0f;
}

// 提取频域特征
void extract_freq_domain_features(float *fft_mag, uint16_t n, FreqDomainFeatures *features) {
float sum_amp = 0.0f, sum_freq_amp = 0.0f;
float max_amp = 0.0f;
uint16_t max_idx = 0;

for (uint16_t i = 0; i < n/2; i++) {
sum_amp += fft_mag[i];
sum_freq_amp += i * fft_mag[i];

if (fft_mag[i] > max_amp) {
max_amp = fft_mag[i];
max_idx = i;
}
}

features->peak_freq = (float)max_idx * SAMPLE_RATE / n;
features->peak_amp = max_amp;
features->total_power = sum_amp;
features->centroid = sum_freq_amp / sum_amp * SAMPLE_RATE / n;
}

8.5 完整 STM32 工程整合

将所有算法集成到主程序中，实现1ms 定时采集→环形缓冲→实时处理→串口输出的完整流程。

// main.c
#include "main.h"
#include "i2c.h"
#include "usart.h"
#include "gpio.h"
#include "tim.h"
#include "mpu6050.h"
#include "vibration_dsp.h"
#include <stdio.h>

MPU6050_DataTypeDef mpu_data;
RingBuffer rb;
AccelData window[FFT_POINTS];
float signal_x[FFT_POINTS];
float fft_mag[FFT_POINTS/2];

TimeDomainFeatures td_features;
FreqDomainFeatures fd_features;

// 报警阈值(根据实际情况调整)
#define RMS_THRESHOLD 0.5f
#define KURTOSIS_THRESHOLD 3.0f

void SystemClock_Config(void);

// 定时器中断服务函数(1ms中断一次，用于精确采集)
void TIM2_IRQHandler(void) {
if (TIM_GetITStatus(TIM2, TIM_IT_Update) != RESET) {
TIM_ClearITPendingBit(TIM2, TIM_IT_Update);

// 读取MPU6050数据
if (MPU6050_ReadAccel(&mpu_data) == 0) {
AccelData data = {mpu_data.accel_x, mpu_data.accel_y, mpu_data.accel_z};
RingBuffer_Push(&rb, &data);
}
}
}

int main(void) {
HAL_Init();
SystemClock_Config();
MX_GPIO_Init();
MX_I2C1_Init();
MX_USART1_UART_Init();
MX_TIM2_Init();

// 初始化环形缓冲区
RingBuffer_Init(&rb);

// 初始化MPU6050
if (MPU6050_Init(&hi2c1) != 0) {
printf("MPU6050初始化失败！\\r\\n");
Error_Handler();
} else {
printf("MPU6050初始化成功！\\r\\n");
printf("开始实时振动信号处理…\\r\\n");
}

// 启动定时器中断(1ms)
HAL_TIM_Base_Start_IT(&htim2);

while (1) {
// 检查是否有足够的数据进行处理
if (RingBuffer_ReadWindow(&rb, window, FFT_POINTS) == 0) {
// 提取X轴信号(工业振动通常以径向为主)
for (uint16_t i = 0; i < FFT_POINTS; i++) {
signal_x[i] = window[i].x;
}

// 预处理
detrend(signal_x, FFT_POINTS);
mean_filter(signal_x, FFT_POINTS, 5);

// 提取时域特征
extract_time_domain_features(signal_x, FFT_POINTS, &td_features);

// FFT频谱分析
// 纯C版本: compute_fft_magnitude(signal_x, fft_mag, FFT_POINTS);
// CMSIS-DSP优化版本:
arm_fft_magnitude(signal_x, fft_mag, FFT_POINTS);

// 提取频域特征
extract_freq_domain_features(fft_mag, FFT_POINTS, &fd_features);

// 简单阈值报警
if (td_features.rms > RMS_THRESHOLD || td_features.kurtosis > KURTOSIS_THRESHOLD) {
HAL_GPIO_WritePin(GPIOA, GPIO_PIN_0, GPIO_PIN_SET); // 点亮报警LED
HAL_GPIO_WritePin(GPIOB, GPIO_PIN_0, GPIO_PIN_SET); // 蜂鸣器报警
} else {
HAL_GPIO_WritePin(GPIOA, GPIO_PIN_0, GPIO_PIN_RESET);
HAL_GPIO_WritePin(GPIOB, GPIO_PIN_0, GPIO_PIN_RESET);
}

// 串口输出处理结果(只输出特征值，不输出原始数据)
printf("%.4f,%.4f,%.4f,%.1f,%.4f\\r\\n",
td_features.rms,
td_features.kurtosis,
fd_features.peak_amp,
fd_features.peak_freq,
fd_features.centroid);
}
}
}

// 重定向printf到USART1
int fputc(int ch, FILE *f) {
HAL_UART_Transmit(&huart1, (uint8_t *)&ch, 1, 100);
return ch;
}

8.6 性能测试与优化建议

STM32F103C8T6 性能测试结果

结论：使用 CMSIS-DSP 库后，512 点 FFT 的处理时间从 85ms 降至 12ms，完全满足 1000Hz 采样率的实时处理要求。

进一步优化建议

九、工业现场工程经验总结

这是我在多个工业项目中总结的宝贵经验，也是很多教程不会告诉你的内容：

十、总结与下期预告

本期我们从零开始搭建了一套基于 STM32+MPU6050 的工业振动数据采集系统，深入讲解了振动信号的预处理、滤波降噪和频谱分析方法。特别地，我们提供了完整的 C 语言实现，将所有信号处理算法直接运行在 STM32 嵌入式端，这是工业现场的标准做法。最后我们制作了可用于模型训练的工业振动故障数据集。

通过本期的学习，你应该已经掌握了：

• 工业振动检测的基本原理和应用场景
• MPU6050 传感器的驱动和数据采集方法
• 振动信号的预处理和滤波降噪技术
• FFT 频谱分析的原理和实现（上位机 Python 版 + 嵌入式 C 版）
• 工业振动故障数据集的制作方法
• 嵌入式端信号处理的工程实践和优化技巧

下期预告：第 6 期《工业振动检测：异常识别模型与故障预警部署》

在下一期中，我们将基于本期制作的数据集和嵌入式端提取的特征，训练两种常用的异常检测模型：