企业官网建设流程全解析

简介：这份MATLAB随机共振代码包面向微弱信号检测与非线性动力学方向的学习者，可用于理解随机共振如何借助噪声增强微弱信号的可检测性，并重点演示变尺度随机共振的实现思路。压缩包内共4个M文件，大小约1KB，代码短小但流程完整，涵盖系统动力学建模、噪声生成、变尺度随机共振数值求解以及输出信号的特征分析等关键环节。已有534人浏览学习这份资源。通过运行和研读这些脚本，读者可以掌握MATLAB中高斯白噪声生成、ode45常微分方程求解、fft频域分析和plot结果可视化等常用操作，体会参数变化对随机共振效果的影响；同时可结合输出信噪比分析脚本学习计算信噪比与评价检测性能的方法。对于正在入门随机共振仿真、或希望将相关方法迁移到轴承故障诊断及弱周期信号提取场景的读者，都具有直接参考价值。

1. 随机共振不是消噪：为什么弱信号检测要靠它

噪声越强，输出信噪比反而越高——这是随机共振最反直觉的地方。它并不是把噪声滤掉，而是让双稳系统的阱间跃迁与弱信号周期同步，把噪声能量“借”给信号频点，让淹没在强噪声里的微弱周期信号重新出现在功率谱上。滚动轴承故障诊断、水声信号检测、生物电信号提取这类低信噪比场景，经常要靠它兜底。MATLAB 做这件事成本极低：一个双稳系统微分方程加一段randn高斯噪声，十几行就能复现，谱估计、扫参、画图都有现成函数。但教材代码通常在超低频率下才有效，实测信号是 kHz 量级，必须引入标题里的“变尺度随机”。这篇文章就从经典模型讲到变尺度实现，最后给出一套可以直接套用的调参验证流程。

2. 从朗之万方程到最小可运行脚本：MATLAB 建模与仿真

2.1 双稳系统如何把噪声能量转到信号频率上

随机共振的物理模型常用一维朗之万方程描述

dx/dt = a*x - b*x^3 + A*sin(2*pi*f0*t) + xi(t)

其中a*x - b*x^3是双稳恢复力，对应的势函数为U(x) = -a/2*x^2 + b/4*x^4，两个稳态分别在x = ±sqrt(a/b)，势垒高度ΔU = a^2/(4b)。粒子在势阱内做噪声驱动的往复运动，弱信号则周期性地抬高或压低势垒。当噪声诱导的 Kramers 跃迁率与信号频率匹配时，粒子跳变会锁定在信号相位上，输出频谱中就会出现与信号同频的谱峰。

关键是绝热近似条件：信号频率、信号幅度、噪声强度都要远小于系统参数。绝大多数 MATLAB 教学代码只在f0=0.01Hz量级下演示，这不是巧合，而是理论推导的前提。一旦输入的频率到了几十赫兹以上，经典公式直接失效，这就是后面变尺度方法存在的理由。

输出时域波形很像一个随机方波：粒子长时间停在某个阱底，偶尔在噪声驱动下跳到另一个阱。这种大幅状态切换本质上是系统对微弱周期信号的阶跃响应，不能再用“小幅振荡”的线性思维去理解。

2.2 最小可运行的随机共振 MATLAB 程序

我一般不用ode45而是用欧拉迭代，原因是随机项会让高阶自适应步长算法变得不稳定，而且欧拉法的步长就是采样间隔，物理含义清楚。最小可运行版本如下：

fs = 101; % 采样率，稍高于信号频率即可 T = 100; % 总时长，保证信号周期数足够多 t = 0:1/fs:T; N = length(t); f0 = 0.01; % 信号频率，满足绝热近似的小参数条件 A = 0.3; % 弱信号幅值 a = 1.0; b = 1.0; % 双稳系统参数，势函数 U = -a/2*x^2 + b/4*x^4 D = 0.6; % 噪声强度 x = zeros(size(t)); for n = 1:N-1 s = A * sin(2*pi*f0*t(n)); xi = sqrt(2*D/fs) * randn; % 注意标定系数是 sqrt(2D/fs) x(n+1) = x(n) + (a*x(n) - b*x(n)^3 + s + xi)/fs; end

两个参数容易踩坑。第一，xi必须乘sqrt(2*D/fs)而不是直接乘D。连续白噪声功率谱密度的离散化，幅值标定系数是sqrt(2D/dt)，漏掉这个系数会让实际注入的噪声强度偏离设定值几个数量级，扫描时共振点会完全对不上。第二，欧拉迭代的+xi写在恢复力和信号的同一级，不要把噪声项单独做一次积分，否则等效于给系统额外加了一个随机游走偏置。

运行后看频谱。用pwelch或者直接abs(fft(x)).^2都可以，重点是在f0=0.01Hz附近看到明显谱峰。阱间跳变产生的随机方波会在低频段贡献很强的能量，所以信号频点周围要选足够窄的窗口，避免把跳变噪声算进信号功率。

2.3 用什么指标判断“共振真的发生”

单看时域波形很难判断，必须量化。常用的指标是输出信噪比：信号频点附近一个窄带窗口内的功率，除以同带宽的噪声功率。

function snr = snr_estimate(pxx, f, f0, bw) idx_sig = abs(f - f0) < bw/2; s = sum(pxx(idx_sig)); n = sum(pxx(~idx_sig)); snr = s / n; end

注意这里~idx_sig把所有频点都当作噪声，得到的并非真正的带内信噪比，只是便于扫描参数时做相对比较。实际应用中应该只取信号频点左右各一个带宽的区间作为噪声窗口，例如abs(f - f0) < bw & ~idx_sig。信噪比增益定义为输出信噪比除以输入信噪比，大于 1 才说明随机共振真的在起作用。

3. 扫描噪声与双稳参数：用信噪比曲线确认共振区间

3.1 扫描噪声强度 D：先升后降的峰就是共振证据

随机共振是否存在，最直接的证据就是“输出信噪比随噪声强度先升后降”。噪声太弱时粒子跃迁不出势垒，信号无法调制状态；噪声太强时跃迁完全随机，信号调制作用被淹没；只有中间的噪声强度让 Kramers 率与信号频率处于同一量级，共振才会出现。

扫描程序如下：

D_list = 0.1:0.05:2.0; snr_list = zeros(size(D_list)); bw = 0.005; for k = 1:numel(D_list) D = D_list(k); x = run_sr(fs, T, f0, A, a, b, D); % 调用上一章的迭代函数 [pxx, f] = pwelch(x, hann(N/8), [], N, fs); snr_list(k) = snr_estimate(pxx, f, f0, bw); end plot(D_list, snr_list, 'LineWidth', 1.5); xlabel('噪声强度 D'); ylabel('输出信噪比');

把迭代部分封装成run_sr函数后，扫描逻辑清晰很多。跑完如果曲线出现单峰，就说明系统在这个参数组合下确实存在随机共振；如果曲线单调上升或单调下降，先检查噪声标定系数有没有漏掉sqrt(2D/fs)，再看信号幅度 A 是否满足小参数条件。

峰值对应的 D 就是当前系统参数下的最优噪声强度。实际信号中噪声强度是固定不可调的，所以扫描 D 的意义不在“调噪声”，而在于验证系统参数是否选在可共振区间内。如果峰值出现在一个非常大的 D 处，说明势垒太高，后面的参数扫描要压低 a 或抬高 b。

3.2 扫描双稳参数 a、b：势垒高度才是核心

固定噪声强度，扫描 a 和 b，本质上是在调节势垒高度ΔU = a^2/(4b)。随机共振发生的必要条件之一是势垒高度与噪声强度同量级：势垒太高粒子出不来，太低信号锁不住。常见经验是ΔU取2D到4D之间。

常用参数调整方向如下表：

二维扫描时可以用surf画出信噪比关于(a, b)的等高面，最佳区域通常在一条曲线附近而不是一个孤点。我一般会先固定 b=1，扫描 a 从 0.2 到 2.0，找到峰值后再微调 b。这样每次只动一个变量，排错更容易。

3.3 实测数据先归一化：量级差太大时直接计算必失败

教科书参数 A=0.3、D=0.6、a=1，碰到实际采集的电压信号动辄几百毫伏甚至几伏，直接代入会直接饱和，输出变成纯方波。处理顺序是先去除直流分量，再按照标准差归一化，让信号幅度和噪声强度都变成 O(1) 量级。

x_raw = x_raw - mean(x_raw); x_norm = x_raw / std(x_raw); D_est = mad(x_raw) / 0.6745; % 用中位绝对偏差估计噪声强度 A_est = bandpower(x_norm, fs, [f0-2 f0+2]); % 粗略估计窄带信号功率

用mad而不是std估计噪声强度，是因为微弱信号对标准差的贡献虽然有但不大，而异常尖峰脉冲会严重抬高std。读取数据文件时也常遇到格式问题：很多采集卡存的是 int16 二进制，用fread读回来如果不按有符号数重新映射，负半周会被折叠成正向大值，归一化后信噪比计算全错。老代码里那些看起来繁琐的字节处理步骤，就是为了修这个问题。

4. 高频信号直接共振会失败：变尺度随机共振的两种实现

4.1 经典模型为什么处理不了 kHz 级信号

绝热近似要求信号频率满足f0 << 1，这对 0.01Hz 的合成信号没问题，换成 500Hz 的实测信号就完全失效。直接代入经典欧拉迭代会发现：步长必须压到1e-4以下，系统参数要么大得离谱，要么输出饱和。

变尺度随机共振的核心思想是时间尺度变换。令tau = alpha * t，朗之万方程变为

dx/dtau = (a/alpha)*x - (b/alpha)*x^3 + (A/alpha)*sin(2*pi*f0/alpha*tau)

信号频率被压缩了alpha倍，系统参数被等效拉伸了alpha倍。也就是说，变尺度处理有两种等价路径：要么对输入信号做时间尺度压缩，让频率落进可共振区；要么保持信号不变，把系统参数 a、b 放大alpha倍。两条路在数学上等价，工程实现和误差特征完全不同。

4.2 路径一：对输入序列做时间尺度压缩（重采样）

我先说适合工程实测的路径：先用interp1对信号做重采样，把高频信号压到 0.01Hz 附近，再做经典随机共振，最后把输出时间轴还原。完整流程如下：

fs_orig = 10000; % 实际采样率 f0_orig = 500; % 目标信号频率 f_target = 0.02; % 压缩后目标频率，落在绝热近似区间 alpha = f0_orig / f_target; % 本例为 25000 L = length(x_raw); t_orig = (0:L-1) / fs_orig; t_resample = (0:ceil(L/alpha)-1) / (fs_orig / alpha); x_resample = interp1(t_orig, x_raw, t_resample, 'spline'); fs_resample = fs_orig / alpha; [x_sr, t_sr] = run_sr(fs_resample, T, f_target, A_est, a, b, D_est); f_recover = f_target * alpha; % 还原为原始频率

重采样之前必须先做抗混叠低通滤波。压缩倍数alpha达到几万倍时，输入信号中高于fs_resample/2的高频分量会被直接折叠进低频区，在输出频谱上表现为信号频点附近的伪峰。建议用lowpass(x_raw, fs_resample/2, fs_orig)先滤一刀。

插值方法的选择也会影响结果。linear速度快但会给信号引入高频转折点；spline保导数光滑，共振曲线更干净，代价是计算量略高。压缩后数据点数只有原来的1/alpha，计算开销反而大幅下降。

4.3 路径二：不重采样，直接拉伸系统参数 a、b

如果不方便改采样率，可以用参数拉伸的等价方法。同样是alpha倍，但作用对象从信号变成系统参数：

a_vs = alpha * a0; b_vs = alpha * b0; [x_sr, t] = run_sr(fs_orig, T, f0_orig, A_est, a_vs, b_vs, D_est);

参数拉伸省去了插值误差，也没有重采样带来的混叠风险，但代价是双稳系统刚性显著增强。alpha = 25000时，欧拉法的稳定步长会变得很小，运行速度反而可能比重采样慢。这时我会改用ode15s或者降采样后再拉伸，并严格检查数值振荡。

两种路径的选择可以参考下表：

实际处理时，我习惯先用参数拉伸快速验证目标频率处是否有共振响应，确认后改用重采样路径做精细的谱分析。两步走能避免把插值误差误判成共振失效。

4.4 变尺度后峰值频偏的排查方法

做完变尺度处理，输出谱峰不一定精确落在f0_recover上，非线性系统会引入一定频偏。允许偏差通常取信号带宽的十倍以内。若偏差太大，先检查alpha是否计算错误，再看重采样前低通截止频率是不是压得太狠，最后检查插值方法是否引入了额外调幅。这四步按顺序排查，绝大多数“共振不出来”的问题都能定位到具体环节。

5. 端到端验证与可复现性：如何判断共振真的发生

5.1 指定参数的回归测试

把整个变尺度流程封装好后，我建议先用已知参数做回归测试：设定一个 500Hz 正弦信号，混入强噪声到输入信噪比约 -15dB，跑完变尺度随机共振流程，检查输出信噪比增益是否大于 2，以及谱峰位置是否落在 500Hz 附近。以下代码是测试骨架：

rng(42); % 固定随机种子，保证结果可复现 N = fs_orig * 10; t = (0:N-1) / fs_orig; x_raw = A*sin(2*pi*f0_orig*t) + sqrt(2*D*fs_orig)*randn(1, N);

固定随机种子这一步骤非常关键。随机共振本身依赖噪声，不固定种子的话，每次运行输出信噪比都有波动，参数扫描结果也难复现。正式实验里我会对每个参数点重复 10 次取平均，用误差棒画信噪比曲线。

5.2 三个能自检的小实验

第一个实验是改变 D 看信噪比曲线峰值位置是否移动。固定信号频率和幅度，分别用 D=0.4、0.6、0.8 跑三次扫描，峰值应向 D 增大的方向移动，且移动幅度应与理论 Kramers 率变化一致。如果曲线形状完全不变，大概率是噪声标定有问题。

第二个实验是检查稳态分布。运行结束后画出状态值的直方图：

figure; histogram(x_sr, 100);

双峰位置应接近±sqrt(a/b)。如果只有单峰，说明势垒太高粒子根本没有跃迁过；如果直方图中间堆满，说明噪声太强系统已经退化为单稳。这个“时域状态→概率分布”的检查方式能直观反映系统是否还处于双稳工作区。

第三个实验是状态切换同步性验证。用阈值把输出波形划分成两个状态，统计状态切换时刻与信号周期的相位差。如果状态切换集中在信号上升沿附近，说明随机共振确实发生了周期同步。这种做法本质上是一种二分聚类，也可以用kmeans自动处理，避免人工选阈值带来的主观偏差。

5.3 导出矢量图时注意的细节

验证通过的曲线最终要导出成论文插图。MATLAB 里我一般用下面这组设置：

set(gcf, 'Color', 'w'); set(findall(gcf, '-property', 'FontSize'), 'FontSize', 9); print(gcf, '-depsc2', 'snr_vs_D.eps');

-depsc2生成的是彩色 EPS 矢量图，缩放不会发虚。如果目标期刊要求在线预览用的无损 PNG，用exportgraphics(gcf, 'snr_vs_D.png', 'Resolution', 600)导出；而用ps2eps做 EPS 裁剪时建议关闭虚线平滑选项，否则功率谱上的细峰容易被光栅化掉。

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