简介面向机械故障诊断、微弱信号检测等研究场景的变尺度随机共振MATLAB实现可用于强噪声背景下微弱周期信号的提取与增强尤其适用于轴承故障特征提取、早期故障诊断等工程问题。代码包含变尺度预处理、随机共振求解、系统输出信噪比计算等核心环节并按函数模块组织便于研究者直接调用、替换信号或修改参数进行仿真对比。资源包共6个文件全部为m脚本整体仅3KB结构紧凑适合有一定MATLAB基础且关注随机共振理论应用的读者快速上手。目前已有182人学习下载说明其实用性得到初步验证。通过这套代码使用者可以理解变尺度随机共振的完整处理流程同时获得一套可扩展的算法骨架用于后续科研实验、课程设计或工程验证中的微弱信号处理任务减少重复编写基础代码的时间。1. 随机共振不是滤噪而是把噪声能量搬给信号当待测信号频率较高比如轴承振动时直接用双稳态随机共振做不出效果因为“小参数”条件不满足。于是出现了变尺度随机共振——先把信号按比例压缩到低频区做完随机共振再映射回原频率。标题里的 bcdSR 就是“随机共振 变尺度”的 MATLAB 实现常用在弱信号检测、故障诊断和生物信号处理。下面我用最小可运行代码和参数表格把从双稳态模型到 bcdSR 的完整链路讲清楚同时指出调参中最容易翻车的位置。2. 双稳态随机共振模型与变尺度随机共振的 MATLAB 建模基础2.1 双稳态朗之万方程随机共振的最小物理模型随机共振stochastic resonance并不是把噪声当作干扰去除而是利用噪声驱动一个非线性系统。最经典的非线性系统是双稳态系统势函数为 U(x) -a/2 x^2 b/4 x^4两个势阱位于 x±√(a/b)。对应的朗之万方程为dx/dt a·x - b·x^3 s(t) ξ(t)其中 s(t) 是待检弱信号ξ(t) 是噪声。当信号幅度 A 远小于势阱深度噪声强度又恰好落在某个区间时粒子在两个势阱间的跳跃频率会和信号频率同步信号频率处的功率被放大。这解释了为什么随机共振能在低信噪比下提升输出信噪比。注意 a 和 b 不是任意缩放都能用。理论上的绝热近似要求信号频率远小于系统弛豫速率。用 MATLAB 仿真时离散步长也要小于系统特征时间否则数值解发散。我先用欧拉-丸山法演示最小模型再给出更适合工程的四阶龙格-库塔形式。2.2 变尺度随机共振把“大参数信号”压回小参数区间经典随机共振的有效工作频率通常在 1 kHz 以下实际工程信号如滚动轴承故障特征频率常在 110 kHz直接输入双稳态系统无法激发共振。变尺度随机共振scale-transformed stochastic resonance通过尺度变换将高频信号“压”到双稳态系统能够响应的低频区间处理后再把频率映射回来。常见的实现是“抽取重采样”设定变换因子 R对长度为 N 的原始信号每隔 R 个点抽取一个样本。这样信号频率变为 f0/R等效采样率变为 fs/R。抽取后的信号用双稳态方程求解求解步长取等效采样间隔 hR/fs。最后将压缩域频谱的横轴乘 R就得到恢复后的原始频率轴。这个流程简单、保留噪声结构是 bcdSR 这类代码包的常用内核。2.3 用 MATLAB 先跑一个最小双稳态随机共振程序在碰 bcdSR 之前我习惯先用一段最小程序确认“共振”存在% minimal_sr.m 最小双稳态随机共振演示 fs 10; % 采样率 10 Hz dt 1/fs; % 步长 N 10000; % 点数 t (0:N-1)*dt; f0 0.033; % 信号频率约 1/30 Hz A 0.3; % 信号幅度 noise_std 0.6; % 噪声标准差 x zeros(1, N); s A*sin(2*pi*f0*t) noise_std*randn(1,N); a 1.0; b 1.0; for i 1:N-1 x(i1) x(i) (a*x(i) - b*x(i)^3 s(i))*dt; end figure; semilogy((0:N-1)/N*fs, abs(fft(s))); hold on; semilogy((0:N-1)/N*fs, abs(fft(x))); legend(输入谱,输出谱); xlabel(Hz); ylabel(幅值);代码里 s 由正弦信号和高斯白噪声相加x 是系统输出。为什么噪声标准差取 0.6这个值大约等于势垒高度的一小部分能让粒子在噪声驱动下周期性越阱。运行后注意看输出谱在 0.033 Hz 附近是否比输入谱有明显抬升。如果抬升了说明参数已经落在共振区抬升不明显就先调噪声强度再调 a/b。这个最小例子不涉及变尺度但它给出了后续 bcdSR 的基线行为。变尺度只是在这个基线上加了重采样一步。3. bcdSR 的 MATLAB 实现变尺度流程与四阶龙格库塔代码3.1 bcdSR 的输入输出与处理阶段bcdSR 这个名称在代码包里一般对应 “bistable stochastic resonance with scale transformation”但不同来源的解压包结构差异很大。下面不依赖某个特定源码直接按工程实现的标准三段式来写输入原始信号和采样率先后完成尺度压缩、双稳态求解、频率映射。处理建议不做带通滤波不做基于小波的去噪。因为随机共振依赖噪声能量提前降噪会削弱效果。如果信号有明显直流偏置先减均值如果有趋势项先用 detrend 去趋势再交给随机共振。这里提一个容易忽略的问题抽取重采样会降低输出长度所以输入的原始信号点数不能太少。我一般要求 N 至少是 R 的 200 倍否则抽取后的序列点数不足频谱分辨率和统计稳定性都差。3.2 变尺度重采样与 RK4 求解核心代码以下是可以直接保存成bcdSR_process.m的函数function [y, f_ax] bcdSR_process(sig, fs, R, a, b) % bcdSR_process: 变尺度随机共振处理 % sig : 输入信号 % fs : 原始采样率 % R : 尺度变换因子 % a/b : 双稳态势函数系数 % y : 增强后的输出信号 % f_ax: 输出频谱对应的原始频率轴 N length(sig); M floor(N / R); sig_c sig(1:R:R*M); % 抽取重采样, 频率降为原频率的 1/R h R / fs; % 等效采样间隔 x zeros(1, M); for i 1:M-1 % 四阶龙格-库塔求解双稳态系统 k1 a*x(i) - b*x(i)^3 sig_c(i); k2 a*(x(i)h*k1/2) - b*(x(i)h*k1/2)^3 sig_c(i); k3 a*(x(i)h*k2/2) - b*(x(i)h*k2/2)^3 sig_c(i); k4 a*(x(i)h*k3) - b*(x(i)h*k3)^3 sig_c(i1); x(i1) x(i) h/6 * (k1 2*k2 2*k3 k4); end f_ax (0:M-1) / (M*h) * R; % 压缩域频率乘 R 映射回原频率 y x; end逻辑说明sig_c是抽取后的压缩序列等效采样率是 fs/R。h 是 RK4 的时间步长按R/fs计算保证积分时每步对应一个重采样点。RK4 的 k2、k3 里驱动项仍用sig_c(i)因为 h 本身是重采样后的采样间隔信号在一个采样周期内变化不大如果信号在局部变化很快可以把驱动项换成插值到中间时刻的数值。f_ax的计算是整个函数的重点压缩域频率分辨率是1/(M*h)乘 R 后得到原始频率位置这样输出频谱可以直接和输入信号画在一起。调用示例fs 100; f0 2; R 4; % 原始 2 Hz 信号压缩到 0.5 Hz t (0:19999)/fs; sig 0.3*sin(2*pi*f0*t) randn(size(t)); [y, f_ax] bcdSR_process(sig, fs, R, 1.0, 1.0); plot(f_ax, abs(fft(y)));这个例子用 100 Hz 采样尺度因子 4把 2 Hz 降到 0.5 Hz正好落在双稳态系统较灵敏的区间。y的频率轴在 0100 Hz 内2 Hz 处可以看到共振后的谱峰。3.3 用频域指标评估随机共振效果只凭输出波形变“干净”了并不算数需要量化。我习惯用“局部峰均比”来对比处理前后的信噪比取 f0 处的谱幅值除以 f0 附近一定带宽内的平均谱幅值然后看处理后的值比处理前高了多少。N_in length(sig); S_in abs(fft(sig)); f_in (0:N_in-1)/N_in*fs; S_out abs(fft(y)); % 注意 f_ax 与 y 等长已经是原频率轴 f0 2; df 0.4; % 信号频率和半带宽 i_f0_in find(abs(f_in-f0)min(abs(f_in-f0))); i_f0_out find(abs(f_ax-f0)min(abs(f_ax-f0))); nei_in find(abs(f_in-f0)df); nei_out find(abs(f_ax-f0)df); peak_ratio_in S_in(i_f0_in) / mean(S_in(nei_in)); peak_ratio_out S_out(i_f0_out) / mean(S_out(nei_out)); gain peak_ratio_out / peak_ratio_in; fprintf(峰均比增益%.3f\n, gain);gain大于 1 说明 bcdSR 提升了目标频率的突出程度。注意df不能取得太小否则邻近噪声样本太少也不能取得太大否则把信号能量也平均掉了。比较稳妥的做法是让2*df覆盖 520 根谱线对 100 Hz 采样、频谱分辨率 0.05 Hz 的情况0.20.4 Hz 都能接受。4. 随机共振参数怎么设bcdSR 的调参顺序与三个常见坑4.1 a、b、R 和等效采样间隔的初始设定没有一组参数能通吃所有信号但可以按下面这个经验表定下第一组可运行的数值参数作用典型范围初始建议a势阱深度决定恢复力大小0.151.0b势垒高度决定跳阱难度0.151.0R尺度变换因子控制频率压低倍数2fs/f0使 f0/R 落在 0.020.08 HzhR/fs等效采样间隔由 R 和 fs 导出不手动设噪声强度输入信号中加入的白噪声标准差依信号幅度定先取信号幅度的 25 倍理论上 a1、b1 时系统势垒高度 ΔU a²/(4b) 0.25弛豫时间约 1 秒量级所以对 0.010.1 Hz 的压缩后信号有响应。若你用的是高频原始信号R 要足够大让压缩后的频率落进这个区间。4.2 参数调整的先后顺序和收敛判断正确顺序是先固定 a、b 和 R再扫描噪声强度或输入信号增益找到输出在 f0 处幅值最大的点然后微调 a 或 b必要时再变 R。直接同时改四个参数会陷入假峰很难复现。收敛判断不要只看时域波形是否像正弦要在频谱上确认两个条件目标频率处幅值持续升高且周围噪声底没有整体抬高。如果增益大于 2对大多数工程弱信号检测场景已经够用如果增益小于 1.2说明参数离最优还远继续扫。4.3 坑 1噪声强度不足时随机共振不出现输入信号如果信噪比太高粒子没有足够噪声驱动无法周期性越阱双稳态系统只会像普通线性放大器一样跟随信号输出不会出现共振峰。表现是扫描噪声强度时输出幅值单调递增没有峰值。解决方法是不要对原始信号做超前降噪甚至在输入到 bcdSR 前人为加一点白噪声。随机共振的本质就是“噪声换信噪比”输入完全没有噪声反而无效。工程上可以在一个范围内扫描噪声增益比如 0.5 倍、1 倍、2 倍、4 倍信号幅度找到峰均比最大的一档。4.4 坑 2尺度因子 R 和采样率不匹配产生伪共振很多资料建议直接取R fs / f0把目标频率压缩到 1 Hz。这个值往往不对因为双稳态系统的最灵敏频率未必是 1 Hz。更可靠的做法是先用纯正弦扫频固定 ab1给不同频率的输入观察输出幅值的峰值在哪一段再回头定 R。伪共振的典型特征是增益曲线随 R 变化不连续R 从 4 改到 5 时效果突然消失改回 4 又出现。这不是共振而是抽取使频谱栅栏错位。把原始数据补零做插值或者增加 N可以缓解。4.5 坑 3RK4 步长过大导致数值振荡在 2.3 节的最小程序中用欧拉法可以但在 bcdSR 里我不建议再用欧拉法。双稳态方程在 x±√(a/b) 附近有刚性步长过大时输出会出现高频振荡甚至发散。检查方法很简单看输出里有没有NaN或者相邻点差分的绝对值突然超过信号幅值的 10 倍。如果出现发散优先检查等效采样间隔 hR/fs 是否小于系统弛豫时间。当 a1、b1 时h 应该明显小于 1 秒如果算出来大于 0.2 秒要么减小 a要么降低 R。不要为了提高频率压缩倍数而无限制增大 R否则数值稳定性会先崩掉。5. 用 bcdSR 处理弱信号后的验证技巧频谱对比与自检5.1 频谱对比法把输入输出画在同一张图里随机共振是否生效最直接的办法是把输入和输出的频谱叠在一起。由于 bcdSR 的输出长度比输入短频率轴要对齐输入用(0:N-1)/N*fs输出直接用函数返回的f_ax。观察目标频率 f0 处是否出现明显凸起同时注意 0 频附近是否被直流抬高。如果直流太高输出谱会在低频段整体翘起来掩盖共振峰这时要先减均值再进 bcdSR。5.2 用合成信号做阳性对照真实数据没有标准答案调参调过头很容易自我欺骗。我的做法是构造一个已知答案的合成信号一个正弦波加高斯白噪声频率和幅度都已知。然后把它送进 bcdSR重复 10 次不同随机噪声种子看峰均比增益的中位数。中位数大于 1.5 才继续用真实数据否则参数没调对先回第 4 章扫参。5.3 快速自检脚本% verify_bcdSR.m 快速自检 fs 100; f0 2; R 4; N 20000; t (0:N-1)/fs; rng(1); sig 0.3*sin(2*pi*f0*t) randn(1,N); [y, f_ax] bcdSR_process(sig, fs, R, 1.0, 1.0); S_in abs(fft(sig)); S_out abs(fft(y)); f_in (0:N-1)/N*fs; figure; plot(f_in, S_in, Color, [0.7 0.7 0.7]); hold on; plot(f_ax, S_out, b); xlim([0 5]); legend(输入,bcdSR输出); xlabel(Hz); ylabel(谱幅值);运行后如果看到 2 Hz 处输出谱线的蓝色高度明显超过灰色输入谱说明链路已经通了。最后提醒一点bcdSR 的输出序列比输入短1/R做频谱对照前一定要用返回的f_ax对齐频率否则峰位偏了一个格你会怀疑算法写错了实际只是坐标轴没用对。本文还有配套的精品资源点击获取 
长河君 · 2026/9/14 13:31:28 · 阅读 1.2k 变尺度随机共振的MATLAB实现：从双稳态模型到bcdSR调参指南
简介面向机械故障诊断、微弱信号检测等研究场景的变尺度随机共振MATLAB实现可用于强噪声背景下微弱周期信号的提取与增强尤其适用于轴承故障特征提取、早期故障诊断等工程问题。代码包含变尺度预处理、随机共振求解、系统输出信噪比计算等核心环节并按函数模块组织便于研究者直接调用、替换信号或修改参数进行仿真对比。资源包共6个文件全部为m脚本整体仅3KB结构紧凑适合有一定MATLAB基础且关注随机共振理论应用的读者快速上手。目前已有182人学习下载说明其实用性得到初步验证。通过这套代码使用者可以理解变尺度随机共振的完整处理流程同时获得一套可扩展的算法骨架用于后续科研实验、课程设计或工程验证中的微弱信号处理任务减少重复编写基础代码的时间。1. 随机共振不是滤噪而是把噪声能量搬给信号当待测信号频率较高比如轴承振动时直接用双稳态随机共振做不出效果因为“小参数”条件不满足。于是出现了变尺度随机共振——先把信号按比例压缩到低频区做完随机共振再映射回原频率。标题里的 bcdSR 就是“随机共振 变尺度”的 MATLAB 实现常用在弱信号检测、故障诊断和生物信号处理。下面我用最小可运行代码和参数表格把从双稳态模型到 bcdSR 的完整链路讲清楚同时指出调参中最容易翻车的位置。2. 双稳态随机共振模型与变尺度随机共振的 MATLAB 建模基础2.1 双稳态朗之万方程随机共振的最小物理模型随机共振stochastic resonance并不是把噪声当作干扰去除而是利用噪声驱动一个非线性系统。最经典的非线性系统是双稳态系统势函数为 U(x) -a/2 x^2 b/4 x^4两个势阱位于 x±√(a/b)。对应的朗之万方程为dx/dt a·x - b·x^3 s(t) ξ(t)其中 s(t) 是待检弱信号ξ(t) 是噪声。当信号幅度 A 远小于势阱深度噪声强度又恰好落在某个区间时粒子在两个势阱间的跳跃频率会和信号频率同步信号频率处的功率被放大。这解释了为什么随机共振能在低信噪比下提升输出信噪比。注意 a 和 b 不是任意缩放都能用。理论上的绝热近似要求信号频率远小于系统弛豫速率。用 MATLAB 仿真时离散步长也要小于系统特征时间否则数值解发散。我先用欧拉-丸山法演示最小模型再给出更适合工程的四阶龙格-库塔形式。2.2 变尺度随机共振把“大参数信号”压回小参数区间经典随机共振的有效工作频率通常在 1 kHz 以下实际工程信号如滚动轴承故障特征频率常在 110 kHz直接输入双稳态系统无法激发共振。变尺度随机共振scale-transformed stochastic resonance通过尺度变换将高频信号“压”到双稳态系统能够响应的低频区间处理后再把频率映射回来。常见的实现是“抽取重采样”设定变换因子 R对长度为 N 的原始信号每隔 R 个点抽取一个样本。这样信号频率变为 f0/R等效采样率变为 fs/R。抽取后的信号用双稳态方程求解求解步长取等效采样间隔 hR/fs。最后将压缩域频谱的横轴乘 R就得到恢复后的原始频率轴。这个流程简单、保留噪声结构是 bcdSR 这类代码包的常用内核。2.3 用 MATLAB 先跑一个最小双稳态随机共振程序在碰 bcdSR 之前我习惯先用一段最小程序确认“共振”存在% minimal_sr.m 最小双稳态随机共振演示 fs 10; % 采样率 10 Hz dt 1/fs; % 步长 N 10000; % 点数 t (0:N-1)*dt; f0 0.033; % 信号频率约 1/30 Hz A 0.3; % 信号幅度 noise_std 0.6; % 噪声标准差 x zeros(1, N); s A*sin(2*pi*f0*t) noise_std*randn(1,N); a 1.0; b 1.0; for i 1:N-1 x(i1) x(i) (a*x(i) - b*x(i)^3 s(i))*dt; end figure; semilogy((0:N-1)/N*fs, abs(fft(s))); hold on; semilogy((0:N-1)/N*fs, abs(fft(x))); legend(输入谱,输出谱); xlabel(Hz); ylabel(幅值);代码里 s 由正弦信号和高斯白噪声相加x 是系统输出。为什么噪声标准差取 0.6这个值大约等于势垒高度的一小部分能让粒子在噪声驱动下周期性越阱。运行后注意看输出谱在 0.033 Hz 附近是否比输入谱有明显抬升。如果抬升了说明参数已经落在共振区抬升不明显就先调噪声强度再调 a/b。这个最小例子不涉及变尺度但它给出了后续 bcdSR 的基线行为。变尺度只是在这个基线上加了重采样一步。3. bcdSR 的 MATLAB 实现变尺度流程与四阶龙格库塔代码3.1 bcdSR 的输入输出与处理阶段bcdSR 这个名称在代码包里一般对应 “bistable stochastic resonance with scale transformation”但不同来源的解压包结构差异很大。下面不依赖某个特定源码直接按工程实现的标准三段式来写输入原始信号和采样率先后完成尺度压缩、双稳态求解、频率映射。处理建议不做带通滤波不做基于小波的去噪。因为随机共振依赖噪声能量提前降噪会削弱效果。如果信号有明显直流偏置先减均值如果有趋势项先用 detrend 去趋势再交给随机共振。这里提一个容易忽略的问题抽取重采样会降低输出长度所以输入的原始信号点数不能太少。我一般要求 N 至少是 R 的 200 倍否则抽取后的序列点数不足频谱分辨率和统计稳定性都差。3.2 变尺度重采样与 RK4 求解核心代码以下是可以直接保存成bcdSR_process.m的函数function [y, f_ax] bcdSR_process(sig, fs, R, a, b) % bcdSR_process: 变尺度随机共振处理 % sig : 输入信号 % fs : 原始采样率 % R : 尺度变换因子 % a/b : 双稳态势函数系数 % y : 增强后的输出信号 % f_ax: 输出频谱对应的原始频率轴 N length(sig); M floor(N / R); sig_c sig(1:R:R*M); % 抽取重采样, 频率降为原频率的 1/R h R / fs; % 等效采样间隔 x zeros(1, M); for i 1:M-1 % 四阶龙格-库塔求解双稳态系统 k1 a*x(i) - b*x(i)^3 sig_c(i); k2 a*(x(i)h*k1/2) - b*(x(i)h*k1/2)^3 sig_c(i); k3 a*(x(i)h*k2/2) - b*(x(i)h*k2/2)^3 sig_c(i); k4 a*(x(i)h*k3) - b*(x(i)h*k3)^3 sig_c(i1); x(i1) x(i) h/6 * (k1 2*k2 2*k3 k4); end f_ax (0:M-1) / (M*h) * R; % 压缩域频率乘 R 映射回原频率 y x; end逻辑说明sig_c是抽取后的压缩序列等效采样率是 fs/R。h 是 RK4 的时间步长按R/fs计算保证积分时每步对应一个重采样点。RK4 的 k2、k3 里驱动项仍用sig_c(i)因为 h 本身是重采样后的采样间隔信号在一个采样周期内变化不大如果信号在局部变化很快可以把驱动项换成插值到中间时刻的数值。f_ax的计算是整个函数的重点压缩域频率分辨率是1/(M*h)乘 R 后得到原始频率位置这样输出频谱可以直接和输入信号画在一起。调用示例fs 100; f0 2; R 4; % 原始 2 Hz 信号压缩到 0.5 Hz t (0:19999)/fs; sig 0.3*sin(2*pi*f0*t) randn(size(t)); [y, f_ax] bcdSR_process(sig, fs, R, 1.0, 1.0); plot(f_ax, abs(fft(y)));这个例子用 100 Hz 采样尺度因子 4把 2 Hz 降到 0.5 Hz正好落在双稳态系统较灵敏的区间。y的频率轴在 0100 Hz 内2 Hz 处可以看到共振后的谱峰。3.3 用频域指标评估随机共振效果只凭输出波形变“干净”了并不算数需要量化。我习惯用“局部峰均比”来对比处理前后的信噪比取 f0 处的谱幅值除以 f0 附近一定带宽内的平均谱幅值然后看处理后的值比处理前高了多少。N_in length(sig); S_in abs(fft(sig)); f_in (0:N_in-1)/N_in*fs; S_out abs(fft(y)); % 注意 f_ax 与 y 等长已经是原频率轴 f0 2; df 0.4; % 信号频率和半带宽 i_f0_in find(abs(f_in-f0)min(abs(f_in-f0))); i_f0_out find(abs(f_ax-f0)min(abs(f_ax-f0))); nei_in find(abs(f_in-f0)df); nei_out find(abs(f_ax-f0)df); peak_ratio_in S_in(i_f0_in) / mean(S_in(nei_in)); peak_ratio_out S_out(i_f0_out) / mean(S_out(nei_out)); gain peak_ratio_out / peak_ratio_in; fprintf(峰均比增益%.3f\n, gain);gain大于 1 说明 bcdSR 提升了目标频率的突出程度。注意df不能取得太小否则邻近噪声样本太少也不能取得太大否则把信号能量也平均掉了。比较稳妥的做法是让2*df覆盖 520 根谱线对 100 Hz 采样、频谱分辨率 0.05 Hz 的情况0.20.4 Hz 都能接受。4. 随机共振参数怎么设bcdSR 的调参顺序与三个常见坑4.1 a、b、R 和等效采样间隔的初始设定没有一组参数能通吃所有信号但可以按下面这个经验表定下第一组可运行的数值参数作用典型范围初始建议a势阱深度决定恢复力大小0.151.0b势垒高度决定跳阱难度0.151.0R尺度变换因子控制频率压低倍数2fs/f0使 f0/R 落在 0.020.08 HzhR/fs等效采样间隔由 R 和 fs 导出不手动设噪声强度输入信号中加入的白噪声标准差依信号幅度定先取信号幅度的 25 倍理论上 a1、b1 时系统势垒高度 ΔU a²/(4b) 0.25弛豫时间约 1 秒量级所以对 0.010.1 Hz 的压缩后信号有响应。若你用的是高频原始信号R 要足够大让压缩后的频率落进这个区间。4.2 参数调整的先后顺序和收敛判断正确顺序是先固定 a、b 和 R再扫描噪声强度或输入信号增益找到输出在 f0 处幅值最大的点然后微调 a 或 b必要时再变 R。直接同时改四个参数会陷入假峰很难复现。收敛判断不要只看时域波形是否像正弦要在频谱上确认两个条件目标频率处幅值持续升高且周围噪声底没有整体抬高。如果增益大于 2对大多数工程弱信号检测场景已经够用如果增益小于 1.2说明参数离最优还远继续扫。4.3 坑 1噪声强度不足时随机共振不出现输入信号如果信噪比太高粒子没有足够噪声驱动无法周期性越阱双稳态系统只会像普通线性放大器一样跟随信号输出不会出现共振峰。表现是扫描噪声强度时输出幅值单调递增没有峰值。解决方法是不要对原始信号做超前降噪甚至在输入到 bcdSR 前人为加一点白噪声。随机共振的本质就是“噪声换信噪比”输入完全没有噪声反而无效。工程上可以在一个范围内扫描噪声增益比如 0.5 倍、1 倍、2 倍、4 倍信号幅度找到峰均比最大的一档。4.4 坑 2尺度因子 R 和采样率不匹配产生伪共振很多资料建议直接取R fs / f0把目标频率压缩到 1 Hz。这个值往往不对因为双稳态系统的最灵敏频率未必是 1 Hz。更可靠的做法是先用纯正弦扫频固定 ab1给不同频率的输入观察输出幅值的峰值在哪一段再回头定 R。伪共振的典型特征是增益曲线随 R 变化不连续R 从 4 改到 5 时效果突然消失改回 4 又出现。这不是共振而是抽取使频谱栅栏错位。把原始数据补零做插值或者增加 N可以缓解。4.5 坑 3RK4 步长过大导致数值振荡在 2.3 节的最小程序中用欧拉法可以但在 bcdSR 里我不建议再用欧拉法。双稳态方程在 x±√(a/b) 附近有刚性步长过大时输出会出现高频振荡甚至发散。检查方法很简单看输出里有没有NaN或者相邻点差分的绝对值突然超过信号幅值的 10 倍。如果出现发散优先检查等效采样间隔 hR/fs 是否小于系统弛豫时间。当 a1、b1 时h 应该明显小于 1 秒如果算出来大于 0.2 秒要么减小 a要么降低 R。不要为了提高频率压缩倍数而无限制增大 R否则数值稳定性会先崩掉。5. 用 bcdSR 处理弱信号后的验证技巧频谱对比与自检5.1 频谱对比法把输入输出画在同一张图里随机共振是否生效最直接的办法是把输入和输出的频谱叠在一起。由于 bcdSR 的输出长度比输入短频率轴要对齐输入用(0:N-1)/N*fs输出直接用函数返回的f_ax。观察目标频率 f0 处是否出现明显凸起同时注意 0 频附近是否被直流抬高。如果直流太高输出谱会在低频段整体翘起来掩盖共振峰这时要先减均值再进 bcdSR。5.2 用合成信号做阳性对照真实数据没有标准答案调参调过头很容易自我欺骗。我的做法是构造一个已知答案的合成信号一个正弦波加高斯白噪声频率和幅度都已知。然后把它送进 bcdSR重复 10 次不同随机噪声种子看峰均比增益的中位数。中位数大于 1.5 才继续用真实数据否则参数没调对先回第 4 章扫参。5.3 快速自检脚本% verify_bcdSR.m 快速自检 fs 100; f0 2; R 4; N 20000; t (0:N-1)/fs; rng(1); sig 0.3*sin(2*pi*f0*t) randn(1,N); [y, f_ax] bcdSR_process(sig, fs, R, 1.0, 1.0); S_in abs(fft(sig)); S_out abs(fft(y)); f_in (0:N-1)/N*fs; figure; plot(f_in, S_in, Color, [0.7 0.7 0.7]); hold on; plot(f_ax, S_out, b); xlim([0 5]); legend(输入,bcdSR输出); xlabel(Hz); ylabel(谱幅值);运行后如果看到 2 Hz 处输出谱线的蓝色高度明显超过灰色输入谱说明链路已经通了。最后提醒一点bcdSR 的输出序列比输入短1/R做频谱对照前一定要用返回的f_ax对齐频率否则峰位偏了一个格你会怀疑算法写错了实际只是坐标轴没用对。本文还有配套的精品资源点击获取