简介这是一份面向MATLAB信号处理学习者的微型演示代码包专门解决振动信号由加速度经数值积分转换为速度与位移的核心问题。资源内仅含1个NumIntTest.m脚本以经典梯形积分函数cumtrapz为主线依次示范累积积分、去趋势处理及结果可视化适合正在学习振动分析、结构健康监测或机械故障诊断的初学者快速上手。文件包仅561B代码精简、逻辑清晰便于直接阅读和修改可用于课堂验证或工程预研。已有946人下载学习足见其在教学演示场景中的实用价值。借助该脚本用户可理解不同物理量之间的微分/积分关系掌握对实测加速度信号进行数值积分的完整流程并学会规避低频漂移等常见误差为后续自行处理真实振动数据打下可靠基础。1. 加速度信号积分到底在解决什么问题现场测振动传感器装上去之后拿到的往往是加速度时域波形。设备状态评价却经常要速度和位移转子动平衡看振动速度有效值ISO 10816 的评判区间用的是速度 mm/s轴瓦磨损、低频松动这类故障又得看位移。手里只有加速度信号那就得做积分。问题在于直接对离散加速度序列做累加结果几乎必然漂移——一个微小的直流偏置或者传感器低频噪声经过一次积分变成线性趋势经过两次积分变成抛物线趋势真实振动分量被彻底淹没。这就是振动信号积分处理的典型困境不是“积不出来”而是“积出来的东西没法用”。这篇东西讲两条路时域积分和频域积分把参数怎么设、趋势项怎么去、误差从哪来一次说清楚。适合做设备状态监测、结构测试、数据采集后处理的工程师以及刚把加速度计数据导进 MATLAB 准备做分析的开发人员。2. 时域积分cumtrapz 与趋势项处理的完整流程2.1 为什么不能直接对加速度信号 cumtrapz在 MATLAB 里看到加速度序列第一反应往往是cumtrapz。梯形法累积积分在数值上没有任何问题问题出在积分对象上。加速度传感器输出的信号包含直流偏置、前级放大器的零漂、低频环境干扰。直流偏置哪怕只有 0.001 m/s²一次积分后就会在时间末端形成 0.001×T 的速度偏差时序长度 10 秒偏置就能到 0.01 m/s而一台正常运转的电机振动速度有效值往往也只有 0.52 mm/s。所以时域积分的第一步不是积分本身而是把不该积的分量提前消掉。这个去偏置的过程分为两种层次。第一层是去均值直接a a - mean(a)把整段信号的直流分量抹掉对应频域里把零频分量置零。第二层是去趋势detrend(a)去掉线性趋势但更稳妥的做法是保留一个 DAC 偏移估计的余地先做多项式拟合再减。实际操作里低频漂移往往不是常数而是缓变信号单靠去均值不够需要配合高通滤波把低频段切掉。2.2 一阶积分的最小可跑代码去均值、梯形积分、去趋势function v accel2velocity(a, ts) % a: 加速度时域序列, 单位 m/s^2, 列向量 % ts: 采样间隔, 单位 s % v: 积分后的速度, 单位 m/s a a(:) - mean(a); % 去直流偏置 % 构造时间轴 N length(a); t (0:N-1). * ts; % 梯形法累积积分 v cumtrapz(t, a); % 再次去除积分引入的线性趋势项 v detrend(v, 1); end这段代码的核心逻辑有三步先减均值把模拟链路的直流漂移归零再用cumtrapz(t, a)做数值积分t必须与a等长且从 0 开始最后用detrend(v, 1)把积分过程中由残余低频成分累积出来的线性斜坡拉平。detrend的第二个参数1表示拟合一次多项式并减去等价于线性最小二乘去趋势。如果省略这个参数detrend默认去掉的是线性趋势效果一样。这里有一个值得注意的细节cumtrapz返回的累积积分末点等于全段积分所以速度序列的首个点是 0末点是整段信号的定积分。对于轴承振动这类零均值随机信号定积分本身接近 0趋势项主导了波形形态。去趋势之后v的均值也被拉到 0 附近这符合振动速度的物理含义——纯交变运动没有单向漂移。2.3 二次积分时趋势项怎样被放大的从加速度到位移要做两次积分。如果沿用一阶积分的思路第一层去均值、第二次积分、再去除线性趋势仍会遇到麻烦。原因在于第二次积分会把第一次积分后残留的微小线性趋势变成明显的抛物线形状这已经超出了detrend(v, 1)能处理的范围。function d accel2displacement(a, ts) a a(:) - mean(a); N length(a); t (0:N-1). * ts; v cumtrapz(t, a); v detrend(v, 1); % 第一次积分后去线性趋势 d cumtrapz(t, v); % 二次积分后残余趋势是二阶多项式用二次拟合去除 p polyfit(t, d, 2); d d - polyval(p, t); d d - mean(d); end二次多项式拟合polyfit(t, d, 2)对应的是“位移序列里含有一个常数偏移、一个线性斜坡和一个二次抛物线的叠加”这三项都源自数值积分对残余低频分量的累积不是真实位移。拟合之后做减法就能把波形拽回零均值附近。注意这段代码里没有先做高通滤波所以当地面存在 0.5Hz 以下的强扰动时即便去掉了二次趋势波形中仍然会有大幅低频摆动。此时正确的做法不是继续加大去趋势的阶数——三次、四次多项式拟合会把真实低频振动也一并吃掉——而是转向频域积分用截止频率把不需要的低频段直接掐掉。积分层级趋势项形态消除手段一次积分线性斜坡detrend(v, 1)或一次多项式拟合二次积分二次抛物线二次多项式拟合polyfit(t, d, 2)叠加低频干扰缓慢摆动高通滤波或频域截止去趋势阶数与截止频率要严格区分detrend处理的是“趋势”假设信号是非平稳的缓变过程高通滤波处理的是“频带”假设某一频段以外的能量都是干扰。前者作用于时域阶数一旦超过真实趋势阶数就会损伤波形低频成分后者作用于频域只要有明确的频率边界就不伤真实分量。工程上更推荐后一种做法这就是下一章要展开的频域积分。3. 频域积分用 FFT 把加速度变换成速度和位移的正确写法3.1 积分的频域形式除以 jω 与 (jω)²时域积分是累积求和频域积分则利用傅里叶变换的微分性质——时域里求导对应频域乘以jω反过来积分就对应除以jω。速度是加速度的一次积分所以频率分量A(jω)变成A(jω) / (jω)位移是二次积分对应A(jω) / (jω)²。这个变换的物理含义非常直接低频分量被大幅放大高频分量被压缩所以频域积分天然对低频噪声敏感。也正因如此频域积分必须显式处理零频和接近零频的频点。离散实现时MATLAB 的fft输出是单边排列的双边频谱索引 1 对应 0Hz索引 2 对应fs/N如此递增。要做的就是在构造好频率轴之后把分母jω数组算出来再对 0Hz 和低于截止频率的位置强制置零。分母置零的位置如果直接做除法会出现 NaN 或 Infifft之后整段波形都会变成 NaN。这是频域积分代码里最常见的崩溃点。3.2 频率轴构造、零频处置和低频截止构造频率轴有一处容易写错fft的循环频率从 0 到fs - fs/N共 N 个点。有人习惯用linspace(0, fs, N)但当 N 不是fs的整数倍时linspace的末点会包含fs而第 N 个 FFT 频点对应的是(N-1)*fs/N并不等于 fs这会直接导致频率轴错位。正确写法是f (0:N-1). * fs / N。function [v, d] fft_integrate(a, fs, fc) % a: 加速度时域序列, 单位 m/s^2 % fs: 采样率, 单位 Hz % fc: 低频截止频率, 单位 Hz, 低于该频率的分量置零 % v: 速度, 单位 m/s % d: 位移, 单位 m a a(:) - mean(a); N length(a); Y fft(a); % 双边频谱 f (0:N-1). * fs / N; % 正确的频率轴 omega 2 * pi * f; % 分母构造, 0Hz 处直接赋值为 NaN 占位, 随后整体处理 den_v 1j * omega; den_d (1j * omega).^2; % 低于截止频率的频点和 0Hz 一并置零 zero_mask (f fc) | (f 0); den_v(zero_mask) inf; den_d(zero_mask) inf; V Y ./ den_v; D Y ./ den_d; % 将非法频点的结果强制归零 V(zero_mask) 0; D(zero_mask) 0; v real(ifft(V)); d real(ifft(D)); end这段代码里把所有低于fc的频点以及 0Hz 全部排除在积分之外。分母先设为inf再整体置零比逐点判断更不容易漏。ifft之后取实部是因为数值计算中浮点舍入会在虚部残留极小的数值物理信号只有实部有意义。需要特别强调的是相位。除以jω使每个频率分量相位滞后 90°除以(jω)²使相位滞后 180°。因此频域积分得到的位移与原始加速度是反相关系这在对照解析解或与其他传感器比对时最容易让人困惑——波形明明反了却检查不出代码错误。3.3 fc 取值的边界效应fc不是高通滤波器的陡峭截止而是一个硬置零的边界。实际信号里恰好落在fc附近的频率分量会被突然切断ifft恢复出的时域波形会在首尾出现振铃Gibbs 现象。缓解办法有两个一是把fc取得尽量低于最小关注频率的三分之一二是做一段过渡带让幅度从 0 到 1 线性斜坡过渡而不是硬切。% 过渡带实现: 0.5*fc 到 fc 线性上升 ramp min(1, f / fc * 2); % 0~0.5fc 处为 0, 0.5fc~fc 处线性升至 1 ramp(f 0) 0;将ramp乘到分母数组上替代原来的zero_mask就能把硬截止变成软截止。过渡带越宽时域振铃越轻但代价是真实低频分量也被部分衰减。现场数据里若已知最低故障特征频率是转频的 1 倍频比如 25Hzfc取 25Hz 足够没有必要再留过渡带若分析对象是 0.5Hz 以下的超低频位移fc就不能超过 0.1Hz此时必须保证数据长度足够长让频率分辨率fs/N低于 0.05Hz。4. 从加速度到位移的工程边界参数怎么设、误差从哪里来4.1 采样率、帧长与频率分辨率的关系频域积分的误差来源与参数强相关。频率分辨率由Δf fs / N决定N 是参与单次 FFT 的数据长度。当关注的最低频率是 1Hz而fs2048、N4096 时Δf0.5Hz1Hz 分量落在第二个频点幅度估计本身就有偏差。频率分辨率不够时低频分量的能量会泄漏到相邻频点而这些频点的积分权重1/ω差异很大低频段一个频点的错位就能带来百分之几十的幅值误差。参数典型取法对积分结果的影响fs关注最高频率 × 2.56 以上偏低导致高频混叠污染整个频带N2 的整数次幂或整周期长度决定 ΔfN 过小则 1Hz 以下无法分辨fc最低关注频率 × 0.20.5过高衰减真实分量过低留下漂移窗函数Hanning 或矩形窗矩形窗泄漏大Hanning 幅值略偏小采样率的选择在数据采集阶段就决定了拿到数据后无法补救。分析前先看fs和信号占用的频带来自压电加速度计的数据通常已经内置抗混叠滤波可直接使用来自某些低成本 MEMS 传感器的数据可能需要先看频谱再决定是否重采样。4.2 低频截止频率 fc 的选取原则fc的物理含义是“低于该频率的分量视为无效并丢弃”。它来自两个约束一是传感器自身的低频响应压电型加速度计在 0.5Hz 以下灵敏度急剧下降这里的信号本来就不可信二是数据处理中希望抑制积分漂移低频漂移不管幅值多小经过二次积分后都会变成大幅位移。工程上的取法是先看频谱图。对加速度信号做pwelch观察能量集中区域。旋转机械的振动能量集中在转频及谐频这时fc取转频的 1/10 到 1/5 比较稳。例如转子转速 3000 r/min转频 50Hz则fc设为 5Hz 是合理的。若分析对象是建筑结构或者海浪激励最低关注频率可能低至 0.1Hzfc只能取 0.02Hz 以下此时对数据长度 N 的要求急剧上升。% 用功率谱估计选择 fc 的参考代码 [pxx, f] pwelch(a, hann(8192), 4096, 16384, fs); figure; semilogx(f, 20*log10(pxx)); xlabel(Frequency (Hz)); ylabel(Power (dB)); % 低频段功率不再下降的拐点, 大致就是 fc 可以被抬高的位置选择fc时还要警惕一个误区不要用高通滤波器的截止频率替代fc。高通滤波器比如highpass(a, fc, fs)有过渡带阶数越高过渡带越窄但相位失真越大。频域积分里的fc是硬性的频域置零两者叠加使用时fc要让在高通滤波器的过渡带之外否则截止附近的相位已经被滤波器改变再除以jω得到的波形相位是错的。4.3 泄漏与加窗频域积分最隐蔽的坑对离散信号做 FFT 隐含了周期延拓的假设。如果数据长度不是信号周期的整数倍频谱就会泄漏。泄漏在普通幅值谱里只是让谱峰变胖但到了频域积分这里泄漏出去的旁瓣会被1/ω放大尤其是低频旁瓣放大倍数可以到几十倍。这导致积分结果出现虚假的大幅低频波动。抑制泄漏的手段是加窗。工程上处理振动数据常用 Hanning 窗它对泄漏的抑制效果好但自身也会让波形两端乘以 0 到 1 的包络。加窗后直接积分会造成两端数据失真所以常见做法是对长数据分段加窗每段各自做频域积分然后拼接。拼接时各段起始和结束点被窗函数压低需要重叠 50% 并用交叉淡化合成。% 分段加窗频域积分的骨架, 每段 50% 重叠 seg_len 4096; hop seg_len / 2; win hann(seg_len, periodic); d_out zeros(size(a)); for idx 1:hop:length(a)-seg_len seg a(idx:idxseg_len-1) .* win; [~, d_seg] fft_integrate(seg, fs, fc); % 将 d_seg 累加到输出, 重叠区自然交叉淡化 d_out(idx:idxseg_len-1) d_out(idx:idxseg_len-1) d_seg; end % 以窗函数累加值归一化这里没有做归一化实际使用时需要维护一个win的累加数组再逐点相除。分段处理后帧与帧之间的相位延续性不能保证所以该方法只适合看位移的有效值或包络不适合拼接成连续时域波形用于轴心轨迹分析。5. 用解析解验证积分结果并处理边界与相位校验5.1 正弦激励下的理论对比写完积分函数第一步不是拿真实数据试而是构造一个已知解析解的信号。标准正弦加速度a(t) A sin(2πf₀t)的速度与位移都有精确表达式v(t) -A/(2πf₀) cos(2πf₀t)d(t) -A/(4π²f₀²) sin(2πf₀t)。只要对比积分结果与理论解的有效值误差就能确认代码里频率轴、分母和置零逻辑是否正确。fs 2000; t (0:8191). / fs; f0 20; A 9.81; a A * sin(2*pi*f0*t); [v, d] fft_integrate(a, fs, 1); v_theo -A/(2*pi*f0) * cos(2*pi*f0*t); d_theo -A/(4*pi^2*f0^2) * sin(2*pi*f0*t); rms_v_err rms(v - v_theo) / rms(v_theo); rms_d_err rms(d - d_theo) / rms(d_theo); fprintf(速度RMS误差 %.6f%%, 位移RMS误差 %.6f%%\n, rms_v_err*100, rms_d_err*100);误差在 0.1% 以下说明频域积分核心逻辑正确。若误差偏大优先检查fc是否高于f0——fc若为 20Hz 而信号是 20Hz置零会把整个信号干掉。其次检查f频率轴末点是否达到了fs频率轴多一个点会让相位逐点偏移但幅值误差不会太大。5.2 检查积分结果的三个习惯真实数据没有解析解可以比对就要靠工程一致性判断。习惯一对比速度有效值。用手持测振仪或另外的速度传感器同时采集比较积分得到的速度 RMS 与参考值的差异偏差在 10% 以内算合理。习惯二检查位移波形形态。正常旋转机械的位移时域波形应该是近似正弦且均值接近零若出现明显抛物线包络说明低频截止不够低或者还有趋势项残留。习惯三做带通滤波后再比较。在关注频带内比如 10Hz1000Hz对原始加速度做带通滤波再分别做时域与频域积分两个结果应大致重合。提示真实数据的单位最容易出错。压电加速度计输出常以 mV/g 标定数据采集软件转出来的可能是 V、mV、g、m/s² 四种不同单位。积分前统一换算为 m/s²位移结果才是米否则一条曲线凭空差出 980 倍。遇到.rar分享包里的数据文件先查记录或matfile读出的变量属性确认单位是 g 还是 m/s²再开始整段处理。相位校验可以算一次再验证对正弦加速度积分后把速度波形与理论余弦对比看零点位置是否差 90°。频域积分得到的相位滞后精确等于 90°如果仿真结果偏差超过 5°说明频率轴或分母构造有误不要急着调整滤波参数先回到fft的频率轴定义上找问题。本文还有配套的精品资源点击获取
简介这是一份面向MATLAB信号处理学习者的微型演示代码包专门解决振动信号由加速度经数值积分转换为速度与位移的核心问题。资源内仅含1个NumIntTest.m脚本以经典梯形积分函数cumtrapz为主线依次示范累积积分、去趋势处理及结果可视化适合正在学习振动分析、结构健康监测或机械故障诊断的初学者快速上手。文件包仅561B代码精简、逻辑清晰便于直接阅读和修改可用于课堂验证或工程预研。已有946人下载学习足见其在教学演示场景中的实用价值。借助该脚本用户可理解不同物理量之间的微分/积分关系掌握对实测加速度信号进行数值积分的完整流程并学会规避低频漂移等常见误差为后续自行处理真实振动数据打下可靠基础。1. 加速度信号积分到底在解决什么问题现场测振动传感器装上去之后拿到的往往是加速度时域波形。设备状态评价却经常要速度和位移转子动平衡看振动速度有效值ISO 10816 的评判区间用的是速度 mm/s轴瓦磨损、低频松动这类故障又得看位移。手里只有加速度信号那就得做积分。问题在于直接对离散加速度序列做累加结果几乎必然漂移——一个微小的直流偏置或者传感器低频噪声经过一次积分变成线性趋势经过两次积分变成抛物线趋势真实振动分量被彻底淹没。这就是振动信号积分处理的典型困境不是“积不出来”而是“积出来的东西没法用”。这篇东西讲两条路时域积分和频域积分把参数怎么设、趋势项怎么去、误差从哪来一次说清楚。适合做设备状态监测、结构测试、数据采集后处理的工程师以及刚把加速度计数据导进 MATLAB 准备做分析的开发人员。2. 时域积分cumtrapz 与趋势项处理的完整流程2.1 为什么不能直接对加速度信号 cumtrapz在 MATLAB 里看到加速度序列第一反应往往是cumtrapz。梯形法累积积分在数值上没有任何问题问题出在积分对象上。加速度传感器输出的信号包含直流偏置、前级放大器的零漂、低频环境干扰。直流偏置哪怕只有 0.001 m/s²一次积分后就会在时间末端形成 0.001×T 的速度偏差时序长度 10 秒偏置就能到 0.01 m/s而一台正常运转的电机振动速度有效值往往也只有 0.52 mm/s。所以时域积分的第一步不是积分本身而是把不该积的分量提前消掉。这个去偏置的过程分为两种层次。第一层是去均值直接a a - mean(a)把整段信号的直流分量抹掉对应频域里把零频分量置零。第二层是去趋势detrend(a)去掉线性趋势但更稳妥的做法是保留一个 DAC 偏移估计的余地先做多项式拟合再减。实际操作里低频漂移往往不是常数而是缓变信号单靠去均值不够需要配合高通滤波把低频段切掉。2.2 一阶积分的最小可跑代码去均值、梯形积分、去趋势function v accel2velocity(a, ts) % a: 加速度时域序列, 单位 m/s^2, 列向量 % ts: 采样间隔, 单位 s % v: 积分后的速度, 单位 m/s a a(:) - mean(a); % 去直流偏置 % 构造时间轴 N length(a); t (0:N-1). * ts; % 梯形法累积积分 v cumtrapz(t, a); % 再次去除积分引入的线性趋势项 v detrend(v, 1); end这段代码的核心逻辑有三步先减均值把模拟链路的直流漂移归零再用cumtrapz(t, a)做数值积分t必须与a等长且从 0 开始最后用detrend(v, 1)把积分过程中由残余低频成分累积出来的线性斜坡拉平。detrend的第二个参数1表示拟合一次多项式并减去等价于线性最小二乘去趋势。如果省略这个参数detrend默认去掉的是线性趋势效果一样。这里有一个值得注意的细节cumtrapz返回的累积积分末点等于全段积分所以速度序列的首个点是 0末点是整段信号的定积分。对于轴承振动这类零均值随机信号定积分本身接近 0趋势项主导了波形形态。去趋势之后v的均值也被拉到 0 附近这符合振动速度的物理含义——纯交变运动没有单向漂移。2.3 二次积分时趋势项怎样被放大的从加速度到位移要做两次积分。如果沿用一阶积分的思路第一层去均值、第二次积分、再去除线性趋势仍会遇到麻烦。原因在于第二次积分会把第一次积分后残留的微小线性趋势变成明显的抛物线形状这已经超出了detrend(v, 1)能处理的范围。function d accel2displacement(a, ts) a a(:) - mean(a); N length(a); t (0:N-1). * ts; v cumtrapz(t, a); v detrend(v, 1); % 第一次积分后去线性趋势 d cumtrapz(t, v); % 二次积分后残余趋势是二阶多项式用二次拟合去除 p polyfit(t, d, 2); d d - polyval(p, t); d d - mean(d); end二次多项式拟合polyfit(t, d, 2)对应的是“位移序列里含有一个常数偏移、一个线性斜坡和一个二次抛物线的叠加”这三项都源自数值积分对残余低频分量的累积不是真实位移。拟合之后做减法就能把波形拽回零均值附近。注意这段代码里没有先做高通滤波所以当地面存在 0.5Hz 以下的强扰动时即便去掉了二次趋势波形中仍然会有大幅低频摆动。此时正确的做法不是继续加大去趋势的阶数——三次、四次多项式拟合会把真实低频振动也一并吃掉——而是转向频域积分用截止频率把不需要的低频段直接掐掉。积分层级趋势项形态消除手段一次积分线性斜坡detrend(v, 1)或一次多项式拟合二次积分二次抛物线二次多项式拟合polyfit(t, d, 2)叠加低频干扰缓慢摆动高通滤波或频域截止去趋势阶数与截止频率要严格区分detrend处理的是“趋势”假设信号是非平稳的缓变过程高通滤波处理的是“频带”假设某一频段以外的能量都是干扰。前者作用于时域阶数一旦超过真实趋势阶数就会损伤波形低频成分后者作用于频域只要有明确的频率边界就不伤真实分量。工程上更推荐后一种做法这就是下一章要展开的频域积分。3. 频域积分用 FFT 把加速度变换成速度和位移的正确写法3.1 积分的频域形式除以 jω 与 (jω)²时域积分是累积求和频域积分则利用傅里叶变换的微分性质——时域里求导对应频域乘以jω反过来积分就对应除以jω。速度是加速度的一次积分所以频率分量A(jω)变成A(jω) / (jω)位移是二次积分对应A(jω) / (jω)²。这个变换的物理含义非常直接低频分量被大幅放大高频分量被压缩所以频域积分天然对低频噪声敏感。也正因如此频域积分必须显式处理零频和接近零频的频点。离散实现时MATLAB 的fft输出是单边排列的双边频谱索引 1 对应 0Hz索引 2 对应fs/N如此递增。要做的就是在构造好频率轴之后把分母jω数组算出来再对 0Hz 和低于截止频率的位置强制置零。分母置零的位置如果直接做除法会出现 NaN 或 Infifft之后整段波形都会变成 NaN。这是频域积分代码里最常见的崩溃点。3.2 频率轴构造、零频处置和低频截止构造频率轴有一处容易写错fft的循环频率从 0 到fs - fs/N共 N 个点。有人习惯用linspace(0, fs, N)但当 N 不是fs的整数倍时linspace的末点会包含fs而第 N 个 FFT 频点对应的是(N-1)*fs/N并不等于 fs这会直接导致频率轴错位。正确写法是f (0:N-1). * fs / N。function [v, d] fft_integrate(a, fs, fc) % a: 加速度时域序列, 单位 m/s^2 % fs: 采样率, 单位 Hz % fc: 低频截止频率, 单位 Hz, 低于该频率的分量置零 % v: 速度, 单位 m/s % d: 位移, 单位 m a a(:) - mean(a); N length(a); Y fft(a); % 双边频谱 f (0:N-1). * fs / N; % 正确的频率轴 omega 2 * pi * f; % 分母构造, 0Hz 处直接赋值为 NaN 占位, 随后整体处理 den_v 1j * omega; den_d (1j * omega).^2; % 低于截止频率的频点和 0Hz 一并置零 zero_mask (f fc) | (f 0); den_v(zero_mask) inf; den_d(zero_mask) inf; V Y ./ den_v; D Y ./ den_d; % 将非法频点的结果强制归零 V(zero_mask) 0; D(zero_mask) 0; v real(ifft(V)); d real(ifft(D)); end这段代码里把所有低于fc的频点以及 0Hz 全部排除在积分之外。分母先设为inf再整体置零比逐点判断更不容易漏。ifft之后取实部是因为数值计算中浮点舍入会在虚部残留极小的数值物理信号只有实部有意义。需要特别强调的是相位。除以jω使每个频率分量相位滞后 90°除以(jω)²使相位滞后 180°。因此频域积分得到的位移与原始加速度是反相关系这在对照解析解或与其他传感器比对时最容易让人困惑——波形明明反了却检查不出代码错误。3.3 fc 取值的边界效应fc不是高通滤波器的陡峭截止而是一个硬置零的边界。实际信号里恰好落在fc附近的频率分量会被突然切断ifft恢复出的时域波形会在首尾出现振铃Gibbs 现象。缓解办法有两个一是把fc取得尽量低于最小关注频率的三分之一二是做一段过渡带让幅度从 0 到 1 线性斜坡过渡而不是硬切。% 过渡带实现: 0.5*fc 到 fc 线性上升 ramp min(1, f / fc * 2); % 0~0.5fc 处为 0, 0.5fc~fc 处线性升至 1 ramp(f 0) 0;将ramp乘到分母数组上替代原来的zero_mask就能把硬截止变成软截止。过渡带越宽时域振铃越轻但代价是真实低频分量也被部分衰减。现场数据里若已知最低故障特征频率是转频的 1 倍频比如 25Hzfc取 25Hz 足够没有必要再留过渡带若分析对象是 0.5Hz 以下的超低频位移fc就不能超过 0.1Hz此时必须保证数据长度足够长让频率分辨率fs/N低于 0.05Hz。4. 从加速度到位移的工程边界参数怎么设、误差从哪里来4.1 采样率、帧长与频率分辨率的关系频域积分的误差来源与参数强相关。频率分辨率由Δf fs / N决定N 是参与单次 FFT 的数据长度。当关注的最低频率是 1Hz而fs2048、N4096 时Δf0.5Hz1Hz 分量落在第二个频点幅度估计本身就有偏差。频率分辨率不够时低频分量的能量会泄漏到相邻频点而这些频点的积分权重1/ω差异很大低频段一个频点的错位就能带来百分之几十的幅值误差。参数典型取法对积分结果的影响fs关注最高频率 × 2.56 以上偏低导致高频混叠污染整个频带N2 的整数次幂或整周期长度决定 ΔfN 过小则 1Hz 以下无法分辨fc最低关注频率 × 0.20.5过高衰减真实分量过低留下漂移窗函数Hanning 或矩形窗矩形窗泄漏大Hanning 幅值略偏小采样率的选择在数据采集阶段就决定了拿到数据后无法补救。分析前先看fs和信号占用的频带来自压电加速度计的数据通常已经内置抗混叠滤波可直接使用来自某些低成本 MEMS 传感器的数据可能需要先看频谱再决定是否重采样。4.2 低频截止频率 fc 的选取原则fc的物理含义是“低于该频率的分量视为无效并丢弃”。它来自两个约束一是传感器自身的低频响应压电型加速度计在 0.5Hz 以下灵敏度急剧下降这里的信号本来就不可信二是数据处理中希望抑制积分漂移低频漂移不管幅值多小经过二次积分后都会变成大幅位移。工程上的取法是先看频谱图。对加速度信号做pwelch观察能量集中区域。旋转机械的振动能量集中在转频及谐频这时fc取转频的 1/10 到 1/5 比较稳。例如转子转速 3000 r/min转频 50Hz则fc设为 5Hz 是合理的。若分析对象是建筑结构或者海浪激励最低关注频率可能低至 0.1Hzfc只能取 0.02Hz 以下此时对数据长度 N 的要求急剧上升。% 用功率谱估计选择 fc 的参考代码 [pxx, f] pwelch(a, hann(8192), 4096, 16384, fs); figure; semilogx(f, 20*log10(pxx)); xlabel(Frequency (Hz)); ylabel(Power (dB)); % 低频段功率不再下降的拐点, 大致就是 fc 可以被抬高的位置选择fc时还要警惕一个误区不要用高通滤波器的截止频率替代fc。高通滤波器比如highpass(a, fc, fs)有过渡带阶数越高过渡带越窄但相位失真越大。频域积分里的fc是硬性的频域置零两者叠加使用时fc要让在高通滤波器的过渡带之外否则截止附近的相位已经被滤波器改变再除以jω得到的波形相位是错的。4.3 泄漏与加窗频域积分最隐蔽的坑对离散信号做 FFT 隐含了周期延拓的假设。如果数据长度不是信号周期的整数倍频谱就会泄漏。泄漏在普通幅值谱里只是让谱峰变胖但到了频域积分这里泄漏出去的旁瓣会被1/ω放大尤其是低频旁瓣放大倍数可以到几十倍。这导致积分结果出现虚假的大幅低频波动。抑制泄漏的手段是加窗。工程上处理振动数据常用 Hanning 窗它对泄漏的抑制效果好但自身也会让波形两端乘以 0 到 1 的包络。加窗后直接积分会造成两端数据失真所以常见做法是对长数据分段加窗每段各自做频域积分然后拼接。拼接时各段起始和结束点被窗函数压低需要重叠 50% 并用交叉淡化合成。% 分段加窗频域积分的骨架, 每段 50% 重叠 seg_len 4096; hop seg_len / 2; win hann(seg_len, periodic); d_out zeros(size(a)); for idx 1:hop:length(a)-seg_len seg a(idx:idxseg_len-1) .* win; [~, d_seg] fft_integrate(seg, fs, fc); % 将 d_seg 累加到输出, 重叠区自然交叉淡化 d_out(idx:idxseg_len-1) d_out(idx:idxseg_len-1) d_seg; end % 以窗函数累加值归一化这里没有做归一化实际使用时需要维护一个win的累加数组再逐点相除。分段处理后帧与帧之间的相位延续性不能保证所以该方法只适合看位移的有效值或包络不适合拼接成连续时域波形用于轴心轨迹分析。5. 用解析解验证积分结果并处理边界与相位校验5.1 正弦激励下的理论对比写完积分函数第一步不是拿真实数据试而是构造一个已知解析解的信号。标准正弦加速度a(t) A sin(2πf₀t)的速度与位移都有精确表达式v(t) -A/(2πf₀) cos(2πf₀t)d(t) -A/(4π²f₀²) sin(2πf₀t)。只要对比积分结果与理论解的有效值误差就能确认代码里频率轴、分母和置零逻辑是否正确。fs 2000; t (0:8191). / fs; f0 20; A 9.81; a A * sin(2*pi*f0*t); [v, d] fft_integrate(a, fs, 1); v_theo -A/(2*pi*f0) * cos(2*pi*f0*t); d_theo -A/(4*pi^2*f0^2) * sin(2*pi*f0*t); rms_v_err rms(v - v_theo) / rms(v_theo); rms_d_err rms(d - d_theo) / rms(d_theo); fprintf(速度RMS误差 %.6f%%, 位移RMS误差 %.6f%%\n, rms_v_err*100, rms_d_err*100);误差在 0.1% 以下说明频域积分核心逻辑正确。若误差偏大优先检查fc是否高于f0——fc若为 20Hz 而信号是 20Hz置零会把整个信号干掉。其次检查f频率轴末点是否达到了fs频率轴多一个点会让相位逐点偏移但幅值误差不会太大。5.2 检查积分结果的三个习惯真实数据没有解析解可以比对就要靠工程一致性判断。习惯一对比速度有效值。用手持测振仪或另外的速度传感器同时采集比较积分得到的速度 RMS 与参考值的差异偏差在 10% 以内算合理。习惯二检查位移波形形态。正常旋转机械的位移时域波形应该是近似正弦且均值接近零若出现明显抛物线包络说明低频截止不够低或者还有趋势项残留。习惯三做带通滤波后再比较。在关注频带内比如 10Hz1000Hz对原始加速度做带通滤波再分别做时域与频域积分两个结果应大致重合。提示真实数据的单位最容易出错。压电加速度计输出常以 mV/g 标定数据采集软件转出来的可能是 V、mV、g、m/s² 四种不同单位。积分前统一换算为 m/s²位移结果才是米否则一条曲线凭空差出 980 倍。遇到.rar分享包里的数据文件先查记录或matfile读出的变量属性确认单位是 g 还是 m/s²再开始整段处理。相位校验可以算一次再验证对正弦加速度积分后把速度波形与理论余弦对比看零点位置是否差 90°。频域积分得到的相位滞后精确等于 90°如果仿真结果偏差超过 5°说明频率轴或分母构造有误不要急着调整滤波参数先回到fft的频率轴定义上找问题。本文还有配套的精品资源点击获取
