MATLAB构建滚动轴承与齿轮复合故障仿真信号：从原理到工程实践

简介：本资源是一套面向机械故障诊断研究者与信号处理初学者的MATLAB复合故障仿真工具包，聚焦滚动轴承与齿轮两类关键部件同时发生故障的建模与信号生成问题，有效支撑故障机理分析、诊断算法验证等科研与工程实践。压缩包共7个文件，含3个核心MATLAB脚本（Compound_fault_simulation_signal.m、Envelope.m、PinPu.m）用于故障信号合成与包络解调，2幅JPG图（时域波形、包络谱）和1张PNG公式图直观呈现仿真结果与数学模型，另附1篇CAJ中文核心论文供方法参考；整体大小3.71MB，轻量易用。已有2262人学习下载，用户可直接运行脚本复现复合振动信号，结合包络谱分析提取双重故障特征频率，并借助论文思路延伸至辛几何模态分解等先进方法应用，是理解多源故障耦合机制与开展算法对比实验的实用起点。

1. 项目概述：为什么我们需要复合故障仿真信号？

在设备状态监测与故障诊断领域，滚动轴承和齿轮是旋转机械中最核心、也最易损的部件。过去，我们研究单一故障（比如一个轴承外圈出现点蚀，或者一个齿轮断了一个齿）的振动信号特征，已经形成了一套相对成熟的分析方法。但在真实的工业现场，情况往往复杂得多。一台长期运行的减速机或齿轮箱，轴承和齿轮长期处于耦合工作状态，磨损和疲劳是同步发生的。很可能在轴承出现早期剥落的同时，齿轮也因为齿面磨损产生了调制现象。这时采集到的振动信号，是多种故障源信号的混合体，我们称之为“复合故障信号”。

如果我们的诊断算法只针对单一故障进行训练和测试，一旦遇到这种复合故障场景，其识别准确率往往会大幅下降，甚至产生误判。这就好比一个只学过识别“猫”和“狗”的AI，突然看到一张“猫和狗在一起”的图片，它可能就懵了。因此，为了开发更鲁棒、更贴近实际的智能诊断算法，我们必须有能力在实验室环境下，高保真地模拟出滚动轴承与齿轮同时发生故障的振动信号。这就是本项目——用MATLAB构建复合故障仿真信号程序——的核心价值所在。

这个程序不是一个简单的信号叠加玩具。它需要你深入理解轴承和齿轮各自的故障振动机理，掌握如何用数学模型精确描述这些机理，并最终在MATLAB环境中将它们合成为一个在时域、频域上都与真实信号高度近似的仿真数据。这对于从事故障预测与健康管理（PHM）、机械信号处理、以及人工智能故障诊断算法研究的工程师和学者来说，是一项非常基础且关键的技能。接下来，我将拆解整个构建过程，从原理到代码，分享我踩过的坑和总结的技巧。

2. 核心思路与数学模型构建

仿真信号的逼真度，完全取决于数学模型的准确性。我们不能简单地用几个正弦波加噪声来糊弄，必须从物理机理出发。

2.1 滚动轴承故障振动模型解析

滚动轴承（如深沟球轴承）的局部故障（点蚀、剥落）会产生周期性的冲击振动。其经典模型通常包含以下几个部分：

1. 周期性冲击序列：这是故障特征的核心。假设故障点周期性撞击滚道，会产生一系列衰减振荡。其数学模型可以表示为：y_bearing(t) = sum( A_i * exp(-β*(t - i*T_d)) * sin(2*π*f_n*(t - i*T_d)) )， 其中t - i*T_d > 0。

   • A_i：第i次冲击的幅值，可以设为常数，也可以引入随机波动模拟实际情况。
   • β：衰减系数，决定了冲击振荡衰减的快慢，与系统阻尼有关。
   • f_n：轴承系统的固有频率，这通常不是故障特征频率，而是冲击激发起的结构共振频率，一般在kHz量级。
   • T_d：故障特征周期。这是关键参数，由轴承几何参数和转速决定。例如，外圈故障特征频率（BPFO）的计算公式为：BPFO = (N/2) * f_r * (1 - d/D * cosφ)。其中N是滚珠数，f_r是轴转频，d是滚珠直径，D是节圆直径，φ是接触角。T_d = 1/BPFO。

2. 幅值调制效应：由于故障点位置相对于负载区的变化，冲击的幅值会随着轴的旋转而周期性变化。例如，外圈故障点如果固定在某个位置，当滚珠通过负载区时冲击强，通过非负载区时冲击弱。这会产生一个以轴转频f_r为周期的幅值调制。我们需要用一个低频的周期函数（如(1 + cos(2*π*f_r*t))）去乘以上述冲击序列。

3. 随机滑动效应：在实际中，滚珠存在滑动，故障特征周期会有微小波动。为了更真实，我们可以在每次冲击的时刻i*T_d上引入一个微小的随机时延δ_i，比如δ_i = 0.01*T_d * randn()。

实操心得：初学者最容易犯的错误是混淆f_n和故障特征频率。f_n是高频的共振频（如2000 Hz），它出现在每次冲击的振荡波形里；而故障特征频率（如BPFO=100 Hz）是冲击发生的频率，体现在冲击的间隔上。在频谱上，你会看到在f_n、2f_n等频率周围，出现以BPFO为间隔的边带，这就是“共振解调”技术能提取故障特征的原因。

2.2 齿轮故障振动模型解析

齿轮的振动主要来源于啮合作用，故障（如齿根裂纹、磨损、断齿）会破坏正常的啮合刚度，从而对振动信号产生调制。

1. 正常齿轮啮合振动：可以简化为一个幅值和相位恒定的啮合频率及其谐波组合。y_gear_healthy(t) = sum( A_m * cos(2*π*m*f_m*t + φ_m) )，m=1,2,3...

   • f_m：啮合频率 = 齿轮齿数 * 轴转频 (f_r)。
   • A_m,φ_m：各次谐波的幅值和相位。

2. 故障引入的调制：齿轮故障主要表现为对啮合振动的幅值调制(AM)和频率调制(FM)。

   • 幅值调制（AM）：通常由齿轮的偏心、齿面磨损不均引起，调制频率为故障齿轮所在轴的转频f_r。调制后的信号为：y_gear_AM(t) = (1 + a*cos(2*π*f_r*t)) * y_gear_healthy(t)。
   • 频率调制（FM）：通常由齿距误差、安装误差或裂纹引起，它使啮合频率发生周期性偏移。调制后的信号为：y_gear_FM(t) = A*cos( 2*π*f_m*t + b*sin(2*π*f_r*t) + φ )。其中b是调制指数，决定了频率偏移的幅度。
   • 局部故障（如断齿）：会产生一个以故障齿轮转频f_r为周期的瞬态冲击，这个冲击模型可以与轴承冲击模型类似，但周期更长，振荡频率可能不同。它还会对啮合振动产生一个强烈的、周期性的幅值调制（在断齿经过啮合区时，振动幅值会骤降）。

2.3 复合故障信号的合成与噪声添加

将上述轴承故障信号y_bearing(t)和齿轮故障信号y_gear_faulty(t)线性叠加，就得到了初步的复合故障信号：y_composite(t) = y_bearing(t) + y_gear_faulty(t)。

但这还不够真实。实际采集的信号必然包含噪声。

1. 背景噪声：来自机器其他部件、环境振动等，通常建模为高斯白噪声n_white(t)。
2. 周期性干扰：可能来自电源工频（50/60 Hz）或其谐波，可以用正弦波模拟。
3. 非高斯冲击噪声：模拟随机的外部撞击，可以用泊松过程生成的随机冲击。

因此，最终的仿真信号为：y_final(t) = y_composite(t) + n_white(t) + n_interference(t) + n_impulse(t)。

信噪比（SNR）是一个关键控制参数，用于控制故障特征的明显程度。

3. MATLAB程序实现与核心代码拆解

下面，我们进入实操环节。我将分模块构建这个仿真程序，并附上详细的代码注释和参数说明。

3.1 全局参数设置与基础信号生成

首先，我们需要定义整个仿真系统的基础参数。这部分代码是程序的基石，参数设置不合理，后面全白搭。

%% 1. 全局仿真参数设置 fs = 25600; % 采样频率 (Hz)，根据奈奎斯特定理，应至少为最高分析频率的2倍以上，工程中常取10倍以上。这里25.6kHz足以分析10kHz以下的特征。 T = 2; % 信号总时长 (秒)，时长要能包含足够多的故障周期，一般建议包含几十个周期。 t = 0:1/fs:T-1/fs; % 时间序列 N = length(t); % 采样点数 %% 2. 机械系统参数设置 % 轴系参数 fr = 30; % 轴转频 (Hz)，即转速为1800 RPM。 fn = 4000; % 轴承系统固有频率 (Hz)，这是一个关键参数，需要根据实际轴承-支座结构估算或测试获得。 % 轴承参数 (以6205深沟球轴承为例) N_ball = 9; % 滚珠数量 d = 7.94e-3; % 滚珠直径 (m) D = 39e-3; % 节圆直径 (m) phi = 0; % 接触角 (度)，对于深沟球轴承通常为0。 % 计算轴承外圈故障特征频率 (BPFO) BPFO = (N_ball/2) * fr * (1 - d/D * cosd(phi)); T_bearing = 1/BPFO; % 轴承故障冲击周期 (秒) % 齿轮参数 Z = 25; % 故障齿轮齿数 fm = Z * fr; % 啮合频率 (Hz) T_gear = 1/fr; % 齿轮故障调制周期 (秒)，即轴旋转一周的时间。

注意事项：fn（固有频率）的值非常关键且容易被误设。它并非故障频率，而是冲击激发起的共振频率。如果你不知道具体值，一个经验法则是：对于中小型轴承-支座系统，fn通常在500Hz到5000Hz之间。你可以通过查阅文献或对健康轴承信号进行频谱分析，找到那个在冲击事件中持续出现的共振峰来估计它。

3.2 滚动轴承外圈故障信号生成

我们模拟一个轴承外圈存在剥落故障的情况。

%% 3. 生成滚动轴承外圈故障信号 beta = 800; % 衰减系数，越大衰减越快。需要与fn配合，使冲击振荡在几个周期内衰减到零。 A0 = 5; % 冲击基准幅值 impact_times = 0:T_bearing:T; % 理论冲击发生时刻序列 num_impacts = length(impact_times) - 1; % 冲击次数 % 引入随机滑动：使每次冲击时刻有微小随机偏移，更符合实际 slip_ratio = 0.02; % 滑动率，2%的周期随机波动 actual_impact_times = impact_times(1:num_impacts) + slip_ratio * T_bearing * randn(1, num_impacts); y_bearing = zeros(1, N); for i = 1:num_impacts ti = actual_impact_times(i); % 找到冲击发生时间点对应的采样点索引范围（冲击发生前后一小段时间） idx_range = find(t >= ti & t < ti + 5/beta); % 只计算冲击后5倍衰减时间内的点，提高效率 if ~isempty(idx_range) tau = t(idx_range) - ti; % 生成单次冲击衰减振荡，并加入幅值调制（模拟负载区变化） amplitude_modulation = 0.3 + 0.7 * (1 + cos(2*pi*fr*t(idx_range)))/2; % 调制深度可调 single_impact = A0 * amplitude_modulation .* exp(-beta * tau) .* sin(2*pi*fn * tau); % 将单次冲击添加到总信号中 y_bearing(idx_range) = y_bearing(idx_range) + single_impact; end end % 对轴承信号进行归一化，方便后续调整信噪比 y_bearing = y_bearing / std(y_bearing);

这段代码有几点精妙之处：

1. 效率优化：没有对整个时间序列进行循环卷积，而是只计算每次冲击附近受影响的数据点，这在模拟长信号时能极大提升速度。
2. 幅值调制集成：将幅值调制直接融入到每次冲击的生成循环中，而不是事后再乘，更符合物理过程。
3. 随机滑动：通过actual_impact_times引入了周期波动，避免了理想周期冲击的“栅栏效应”，使仿真频谱更接近真实情况。

3.3 齿轮局部故障（断齿）信号生成

我们模拟一个齿轮存在断齿故障，它会产生周期性冲击并对啮合振动产生强幅值调制。

%% 4. 生成齿轮局部故障（断齿）信号 % 4.1 生成正常的啮合振动（基波和二次谐波） A_m1 = 1.2; A_m2 = 0.6; % 啮合频率基波和二次谐波幅值 phi_m1 = pi/6; phi_m2 = pi/4; % 初始相位 y_gear_healthy = A_m1 * cos(2*pi*fm*t + phi_m1) + A_m2 * cos(2*pi*2*fm*t + phi_m2); % 4.2 模拟断齿引起的周期性冲击（每转一次） % 断齿冲击的模型可能与轴承冲击不同，通常衰减更快，振荡频率可能接近齿轮箱的固有频率。 fn_gear = 1200; % 齿轮箱体固有频率 (Hz)，区别于轴承固有频率。 beta_gear = 1200; % 齿轮冲击衰减系数 A_gear_impact = 8; % 断齿冲击幅值 gear_impact_times = 0:T_gear:T; % 每转一次冲击 num_gear_impacts = length(gear_impact_times) - 1; y_gear_impact = zeros(1, N); for i = 1:num_gear_impacts ti = gear_impact_times(i); idx_range = find(t >= ti & t < ti + 5/beta_gear); if ~isempty(idx_range) tau = t(idx_range) - ti; single_gear_impact = A_gear_impact * exp(-beta_gear * tau) .* sin(2*pi*fn_gear * tau); y_gear_impact(idx_range) = y_gear_impact(idx_range) + single_gear_impact; end end % 4.3 断齿对啮合振动的强幅值调制 % 当断齿经过啮合区时，振动幅值会几乎为零。用一个周期性的矩形窗函数来模拟。 modulation_depth = 0.95; % 调制深度，0.95表示在故障位置振幅降至5%。 % 创建一个周期性的调制函数：故障角位置（假设断齿在0相位角）附近幅值最小。 fault_angle_width = 0.1; % 故障齿对应的角度宽度（归一化到0-1），0.1表示36度角范围内受影响。 angle_series = mod(fr * t, 1); % 将时间转换为轴旋转的周期相位（0-1） % 使用一个平滑的凹陷函数来模拟调制，比矩形窗更真实（可用余弦窗） modulation_func = 1 - modulation_depth * exp(-(angle_series.^2)/(2*(fault_angle_width/3)^2)); % 高斯窗模拟平滑过渡 y_gear_modulated = modulation_func .* y_gear_healthy; % 4.4 合成齿轮故障信号：啮合振动（被调制）+ 断齿冲击 y_gear_fault = y_gear_modulated + y_gear_impact; y_gear_fault = y_gear_fault / std(y_gear_fault); % 归一化

这里的关键在于对“断齿调制”的处理。我没有使用简单的0-1矩形窗，而是用一个平滑的高斯窗函数，这避免了引入不真实的高频成分，使得仿真信号的包络谱分析更合理。

3.4 信号合成、加噪与可视化

现在，我们将轴承故障和齿轮故障信号合二为一，并加入噪声环境。

%% 5. 合成复合故障信号并添加噪声 % 5.1 合成无噪声的复合故障信号 composite_clean = y_bearing + 0.7 * y_gear_fault; % 调整权重，模拟不同故障的严重程度 composite_clean = composite_clean / std(composite_clean); % 再次归一化 % 5.2 添加高斯白噪声（背景噪声） target_snr_db = 0; % 目标信噪比 (dB)，0dB表示信号功率与噪声功率相等。可以调整，负值表示噪声更强。 signal_power = sum(composite_clean.^2) / N; noise_power = signal_power / (10^(target_snr_db/10)); noise_white = sqrt(noise_power) * randn(1, N); % 5.3 添加周期性干扰（例如，电源线干扰） f_interf = 50; % 工频干扰 50 Hz A_interf = 0.3; % 干扰幅值 noise_interference = A_interf * sin(2*pi*f_interf*t); % 5.4 合成最终信号 y_final = composite_clean + noise_white + noise_interference; %% 6. 时域波形可视化 figure('Position', [100, 100, 1200, 800]) subplot(4,2,1) plot(t, y_bearing) title('纯轴承外圈故障信号 (时域)') xlabel('时间 (s)'); ylabel('幅值'); grid on; xlim([0, 0.2]) % 显示前0.2秒 subplot(4,2,2) plot(t, y_gear_fault) title('纯齿轮断齿故障信号 (时域)') xlabel('时间 (s)'); ylabel('幅值'); grid on; xlim([0, 0.2]) subplot(4,2,3) plot(t, composite_clean) title('复合故障信号 (清洁版)') xlabel('时间 (s)'); ylabel('幅值'); grid on; xlim([0, 0.2]) subplot(4,2,4) plot(t, y_final) title('复合故障信号 (最终含噪版)') xlabel('时间 (s)'); ylabel('幅值'); grid on; xlim([0, 0.2]) %% 7. 频域分析（功率谱密度） [psd_bearing, f_b] = pwelch(y_bearing, hanning(4096), 2048, 4096, fs); [psd_gear, f_g] = pwelch(y_gear_fault, hanning(4096), 2048, 4096, fs); [psd_composite, f_c] = pwelch(composite_clean, hanning(4096), 2048, 4096, fs); [psd_final, f_f] = pwelch(y_final, hanning(4096), 2048, 4096, fs); subplot(4,2,5) semilogy(f_b, psd_bearing) title('纯轴承故障信号 (频谱)'); xlabel('频率 (Hz)'); ylabel('功率谱密度'); grid on; xlim([0, fs/2]) subplot(4,2,6) semilogy(f_g, psd_gear) title('纯齿轮故障信号 (频谱)'); xlabel('频率 (Hz)'); ylabel('功率谱密度'); grid on; xlim([0, fs/2]) subplot(4,2,7) semilogy(f_c, psd_composite) title('复合故障信号 (频谱，清洁)'); xlabel('频率 (Hz)'); ylabel('功率谱密度'); grid on; xlim([0, fs/2]) subplot(4,2,8) semilogy(f_f, psd_final) title('复合故障信号 (频谱，含噪)'); xlabel('频率 (Hz)'); ylabel('功率谱密度'); grid on; xlim([0, fs/2])

可视化是验证仿真结果是否合理的必要步骤。通过对比时域波形和频谱图，你可以清晰地看到：

• 轴承故障信号在fn(4000 Hz) 附近有明显的共振带。
• 齿轮故障信号在fm(750 Hz) 和2fm(1500 Hz) 处有清晰的谱线，并且在fr(30 Hz) 间隔处存在边带（调制现象）。
• 复合信号同时包含了以上特征。
• 加入噪声后，低频段的故障特征（如BPFO的谐波）可能被淹没，但高频共振带依然可见，这模拟了实际中常使用高频共振解调技术的原因。

4. 高级分析与特征提取验证

生成了信号只是第一步，我们还需要用标准的诊断方法来验证仿真信号是否包含了预期的故障特征。这是检验模型有效性的关键。

4.1 包络谱分析（用于轴承故障诊断）

包络谱分析是提取轴承周期性冲击特征的最有效工具之一。其原理是先对信号进行带通滤波（通常围绕系统共振频率fn），然后通过希尔伯特变换提取滤波后信号的包络，最后对包络信号进行频谱分析。

%% 8. 对最终复合信号进行包络谱分析（验证轴承故障特征） % 8.1 带通滤波：围绕轴承系统固有频率 fn bp_filt = designfilt('bandpassiir', 'FilterOrder', 8, ... 'HalfPowerFrequency1', fn-500, 'HalfPowerFrequency2', fn+500, ... 'SampleRate', fs); y_filtered = filtfilt(bp_filt, y_final); % 使用零相位滤波filtfilt避免相位失真 % 8.2 希尔伯特变换提取包络 y_envelope = abs(hilbert(y_filtered)); % 8.3 计算包络信号的频谱（即包络谱） [psd_env, f_env] = pwelch(y_envelope, hanning(4096), 2048, 4096, fs); figure; subplot(2,1,1) plot(t, y_envelope); xlim([0, 1]); % 显示包络波形 title('带通滤波后信号的包络'); xlabel('时间 (s)'); ylabel('幅值'); grid on; subplot(2,1,2) plot(f_env, psd_env); xlim([0, 500]); % 重点关注低频段，故障特征频率在此 title('包络谱 (应能观察到轴承故障特征频率BPFO及其谐波)'); xlabel('频率 (Hz)'); ylabel('功率谱密度'); grid on; hold on; % 在图上标记理论BPFO位置 bpfo_line = line([BPFO, BPFO], ylim, 'Color', 'r', 'LineStyle', '--', 'LineWidth', 1.5); text(BPFO*1.05, max(psd_env)*0.9, sprintf('BPFO = %.2f Hz', BPFO), 'Color', 'r'); % 标记其谐波 for k=2:5 line([k*BPFO, k*BPFO], ylim, 'Color', 'r', 'LineStyle', ':', 'LineWidth', 1); end hold off;

如果模型正确，在包络谱的低频部分（通常0-500Hz），你应该能看到在BPFO、2*BPFO、3*BPFO等处出现明显的谱峰。这证明了仿真信号中成功嵌入了轴承外圈故障的周期性冲击特征。

4.2 频谱细化与边带分析（用于齿轮故障诊断）

对于齿轮故障，我们需要在啮合频率fm周围观察是否存在以轴转频fr为间隔的边带，这是幅值/频率调制的典型特征。

%% 9. 频谱细化分析（验证齿轮故障调制边带） % 9.1 计算高分辨率频谱（使用更大的FFT点数） N_fft = 2^nextpow2(10*fs); % 使用10秒数据长度的FFT，提高频率分辨率 [psd_highres, f_highres] = pwelch(y_final, hanning(N_fft/4), N_fft/8, N_fft, fs); % 9.2 聚焦在啮合频率 fm 附近 f_range_low = fm - 5*fr; % 观察fm前后5倍转频的范围 f_range_high = fm + 5*fr; idx_range = find(f_highres >= f_range_low & f_highres <= f_range_high); figure; plot(f_highres(idx_range), 10*log10(psd_highres(idx_range))); % 用dB显示更清晰 title(sprintf('啮合频率fm (%.1f Hz) 附近的细化频谱', fm)); xlabel('频率 (Hz)'); ylabel('功率谱密度 (dB)'); grid on; hold on; % 标记啮合频率fm plot([fm, fm], ylim, 'k--', 'LineWidth', 1.5); text(fm, max(10*log10(psd_highres(idx_range))), ' f_m', 'VerticalAlignment', 'top'); % 标记可能的边带位置 for k = -5:5 if k ~= 0 f_sideband = fm + k * fr; plot([f_sideband, f_sideband], ylim, 'r:', 'LineWidth', 0.8); end end hold off;

在细化频谱图中，如果齿轮存在故障调制，你会在fm谱线的两侧，间隔fr的位置，看到明显的边带谱线。边带的幅值和数量反映了调制的深度和类型。

4.3 时频分析（观察瞬态冲击）

对于断齿这类瞬态冲击事件，时频分析（如短时傅里叶变换STFT）可以直观展示冲击发生的时间及其频率成分。

%% 10. 时频分析（观察齿轮断齿冲击） figure; [s, f_stft, t_stft] = spectrogram(y_final, hanning(512), 256, 512, fs, 'yaxis'); % 使用imagesc绘制时频谱图，更灵活 imagesc(t_stft, f_stft, 10*log10(abs(s))); axis xy; % 确保频率轴方向正确 xlabel('时间 (s)'); ylabel('频率 (Hz)'); title('复合故障信号时频谱图 (STFT)'); colorbar; ylim([0, 3000]); % 聚焦在低频段，观察冲击和调制 % 在图上叠加理论冲击时刻线（齿轮转频） hold on; for i = 0:floor(T/T_gear) plot([i*T_gear, i*T_gear], ylim, 'w--', 'LineWidth', 1); end hold off;

在时频谱图中，你应该能看到垂直方向的“条纹”，其出现间隔约为T_gear(1/30 ≈ 0.033秒)，这对应了齿轮每转一周断齿产生的冲击。同时，在fm频率处应该有一条持续的水平亮线，但其亮度可能会在冲击时刻附近发生变化，这体现了幅值调制。

5. 参数化与批量生成：构建你的仿真数据集

一个优秀的仿真程序不能只是生成一个信号。为了用于机器学习算法训练，我们需要能批量生成不同工况、不同故障程度、不同信噪比的数据。

5.1 将仿真过程封装为函数

我们将核心仿真逻辑封装成一个MATLAB函数，便于调用。

function [signal, params] = generate_composite_fault_signal(params) % GENERATE_COMPOSITE_FAULT_SIGNAL 生成滚动轴承与齿轮复合故障仿真信号 % 输入：params - 结构体，包含所有仿真参数 % 输出：signal - 生成的仿真信号 % params - 返回使用的参数（可能包含默认值） % 设置默认参数 default_params.fs = 25600; default_params.T = 2; default_params.fr = 30; default_params.fn_bearing = 4000; default_params.beta_bearing = 800; default_params.BPFO = 90; % 或通过几何参数计算 default_params.fn_gear = 1200; default_params.beta_gear = 1200; default_params.Z = 25; default_params.fault_severity_bearing = 1.0; % 轴承故障严重程度因子 default_params.fault_severity_gear = 1.0; % 齿轮故障严重程度因子 default_params.SNR_dB = 0; % 信噪比 default_params.add_50Hz_interference = true; % 是否添加50Hz干扰 % 用输入参数覆盖默认值 if nargin < 1 params = struct(); end param_fields = fieldnames(default_params); for i = 1:length(param_fields) field = param_fields{i}; if ~isfield(params, field) params.(field) = default_params.(field); end end % 解包参数（代码略，与前面章节类似，但所有硬编码数字替换为params.xxx） % ... (此处插入前面第3.1至3.4节的代码，但将所有数字改为参数变量) ... % 例如：fs = params.fs; fr = params.fr; 等等。 % 最终返回信号和参数 signal = y_final; end

5.2 设计实验矩阵与批量生成

假设我们要研究不同信噪比和不同故障严重程度组合下的信号特征。

%% 11. 批量生成仿真数据集 % 定义实验变量 snr_levels = [10, 0, -5]; % 信噪比 (dB)，从清晰到嘈杂 bearing_severity_levels = [0.5, 1.0, 2.0]; % 轴承故障严重程度因子 gear_severity_levels = [0.5, 1.0, 2.0]; % 齿轮故障严重程度因子 num_samples_per_case = 50; % 每种工况生成50个样本 sample_length = params.fs * params.T; % 每个样本的长度 % 初始化数据存储（假设用单元数组或结构体数组） dataset = cell(length(snr_levels) * length(bearing_severity_levels) * length(gear_severity_levels) * num_samples_per_case, 1); labels = zeros(size(dataset, 1), 3); % 标签列：[SNR, 轴承严重度, 齿轮严重度] idx = 1; h = waitbar(0, '正在批量生成仿真数据...'); total_iters = length(snr_levels) * length(bearing_severity_levels) * length(gear_severity_levels) * num_samples_per_case; current_iter = 0; for snr = snr_levels for s_b = bearing_severity_levels for s_g = gear_severity_levels for n = 1:num_samples_per_case % 设置当前工况参数 current_params = params; % params是基础参数结构体 current_params.SNR_dB = snr; current_params.fault_severity_bearing = s_b; current_params.fault_severity_gear = s_g; % 可以在每次生成时加入微小的随机参数波动，增加数据多样性 current_params.fr = params.fr * (1 + 0.01*randn()); % 转速微小波动 current_params.BPFO = current_params.BPFO * (1 + 0.01*randn()); % 故障频率波动 % 生成信号 [signal, ~] = generate_composite_fault_signal(current_params); % 存储数据和标签 dataset{idx} = signal; labels(idx, :) = [snr, s_b, s_g]; idx = idx + 1; current_iter = current_iter + 1; waitbar(current_iter/total_iters, h); end end end end close(h); fprintf('批量生成完成！共生成 %d 个样本。\n', length(dataset)); % 可以保存数据集供后续使用 % save('composite_fault_dataset.mat', 'dataset', 'labels', 'params', '-v7.3');

通过这样的批量生成，你就能快速构建一个规模可观、工况覆盖全面的复合故障仿真数据集，用于后续深度学习或机器学习模型的训练和测试。

6. 常见问题、调试技巧与进阶思考

在实际编写和运行这个程序时，你肯定会遇到各种问题。下面是我总结的一些常见坑点和解决思路。

6.1 频谱上找不到预期的故障特征频率？

这是最常见的问题。请按以下步骤排查：

6.2 程序运行速度太慢？

生成长信号或批量生成时，循环可能是瓶颈。优化策略：

• 向量化操作：尽可能用矩阵运算代替循环。例如，冲击序列的生成可以用toeplitz矩阵或卷积实现，但对于变周期和随机滑动，循环更直观。我们的代码已对单次冲击进行了局部计算优化。
• 使用parfor并行循环：在批量生成数据时（第5.2节），最外层的循环可以使用parfor进行并行计算，大幅提升速度。确保你的MATLAB安装了Parallel Computing Toolbox。

  parfor idx = 1:total_samples % 生成信号的代码 end

• 预分配数组：在循环前用zeros()预分配所有数组空间，避免动态增长。

6.3 如何让仿真信号更贴近真实？

1. 引入传递函数：真实的振动信号从故障源传到传感器，会经过复杂的路径，相当于经过一个“传递函数”滤波。你可以在生成故障冲击后，用一个模拟结构传递特性的IIR/FIR滤波器（如低通或带通）对信号进行滤波。
2. 非平稳转速：实际设备转速会有波动。可以将轴转频fr设置为时间的函数fr(t)，然后通过积分得到瞬时相位φ(t) = ∫ 2π fr(t) dt，再用φ(t)替换所有公式中的2π f t。这能模拟升降速过程。
3. 更复杂的噪声模型：除了白噪声，可以加入“有色噪声”（功率谱不平坦）和“非平稳噪声”，模拟实际工况。
4. 多故障点组合：模拟轴承内圈、滚动体故障与齿轮故障的组合，或者多个轴承同时故障。

6.4 这个仿真程序的主要应用场景是什么？

1. 算法开发与验证：为新的故障诊断算法（如深度学习模型、自适应滤波算法、盲源分离算法）提供可控、可重复、带精确标签的测试数据。
2. 特征提取方法研究：研究在复合故障、强噪声背景下，不同特征提取方法（如包络谱、小波变换、循环平稳分析）的有效性和鲁棒性。
3. 教学与理解：直观展示轴承和齿轮故障信号的产生机理、调制边带、共振解调等抽象概念，是理解机械故障诊断原理的绝佳工具。
4. 数据增强：当真实故障数据稀缺时，可以用仿真数据对模型进行预训练或数据增强，提升模型在实际应用中的泛化能力。

构建这样一个复合故障仿真程序，就像为故障诊断研究搭建了一个高保真的“数字风洞”。它允许你在成本极低的情况下，无限次地复现各种复杂的故障工况，从而加速算法的迭代与验证。我强烈建议你在理解上述代码的基础上，动手调整每一个参数，观察信号时域、频域、时频域的变化，这比读任何教科书都更能加深你对机械故障振动特征的理解。当你能够通过调整参数，让仿真信号的频谱图与你手中真实的故障谱图高度相似时，你就真正掌握了这门“手艺”。

本文还有配套的精品资源，点击获取