MATLAB实现，用于模拟和分析归一化随机共振现象。随机共振是一种非线性现象，当弱周期性信号与适当强度的噪声共同作用时，系统输出信噪比会显著提高。
%% 归一化随机共振仿真
% 功能: 实现双稳态系统的归一化随机共振仿真
% 模型: 一维双势阱系统clear; clc; close all;%% 1. 参数设置
% 系统参数
a = 1.0; % 势阱深度参数
b = 1.0; % 势阱宽度参数% 信号参数
A = 0.1; % 信号幅值
f0 = 0.05; % 信号频率 (Hz)
phase = 0; % 信号相位 (弧度)% 噪声参数
D = 0.2; % 噪声强度
noise_type = 'gaussian'; % 噪声类型: 'gaussian', 'uniform', 'colored'% 仿真参数
t_start = 0; % 起始时间
t_end = 1000; % 结束时间
dt = 0.1; % 时间步长
t = t_start:dt:t_end; % 时间向量
N = length(t); % 时间点数% 归一化参数
normalize_output = true; % 是否归一化输出
normalization_method = 'minmax'; % 'minmax' 或 'zscore'%% 2. 生成输入信号
% 周期性信号
signal = A * sin(2*pi*f0*t + phase);% 添加直流偏置（可选）
dc_offset = 0.0;
signal = signal + dc_offset;%% 3. 生成噪声
switch noise_typecase 'gaussian'% 高斯白噪声noise = sqrt(2*D*dt) * randn(1, N);case 'uniform'% 均匀噪声noise = sqrt(12*D*dt) * (rand(1, N) - 0.5);case 'colored'% 有色噪声（一阶低通滤波）raw_noise = sqrt(2*D*dt) * randn(1, N);alpha = 0.9; % 滤波系数colored_noise = zeros(1, N);colored_noise(1) = raw_noise(1);for i = 2:Ncolored_noise(i) = alpha*colored_noise(i-1) + (1-alpha)*raw_noise(i);endnoise = colored_noise;otherwiseerror('未知噪声类型');
end%% 4. 归一化处理
% 归一化输入信号
if normalize_outputswitch normalization_methodcase 'minmax'% 最小-最大归一化signal_norm = (signal - min(signal)) / (max(signal) - min(signal));noise_norm = (noise - min(noise)) / (max(noise) - min(noise));case 'zscore'% Z-score标准化signal_norm = (signal - mean(signal)) / std(signal);noise_norm = (noise - mean(noise)) / std(noise);otherwiseerror('未知归一化方法');end
elsesignal_norm = signal;noise_norm = noise;
end%% 5. 双稳态系统模型
% 初始化状态变量
x = zeros(1, N);
x(1) = 0.1; % 初始条件% 系统动力学方程
% dx/dt = a*x - b*x^3 + signal + noise
for i = 1:N-1dx = (a*x(i) - b*x(i)^3 + signal_norm(i) + noise_norm(i)) * dt;x(i+1) = x(i) + dx;
end%% 6. 归一化输出
if normalize_outputswitch normalization_methodcase 'minmax'output = (x - min(x)) / (max(x) - min(x));case 'zscore'output = (x - mean(x)) / std(x);otherwiseoutput = x;end
elseoutput = x;
end%% 7. 结果可视化
% 时域信号
figure('Name', '时域信号', 'Position', [100, 100, 1200, 800]);
subplot(3,1,1);
plot(t, signal, 'b', 'LineWidth', 1.5);
title('输入信号');
xlabel('时间');
ylabel('幅值');
grid on;subplot(3,1,2);
plot(t, noise, 'r', 'LineWidth', 1.5);
title('噪声');
xlabel('时间');
ylabel('幅值');
grid on;subplot(3,1,3);
plot(t, output, 'k', 'LineWidth', 1.5);
title('系统输出');
xlabel('时间');
ylabel('幅值');
grid on;% 相空间图
figure('Name', '相空间图', 'Position', [100, 100, 800, 600]);
plot(output(1:end-1), diff(output)*dt, '.', 'MarkerSize', 1);
title('相空间轨迹');
xlabel('x');
ylabel('dx/dt');
grid on;% 直方图
figure('Name', '状态分布', 'Position', [100, 100, 800, 600]);
histogram(output, 50, 'Normalization', 'pdf');
hold on;% 理论概率密度函数
x_vals = linspace(min(output), max(output), 100);
pdf_stable1 = (1/sqrt(2*pi)) * exp(-(x_vals - sqrt(a/b)).^2/(2*a/(2*b)));
pdf_stable2 = (1/sqrt(2*pi)) * exp(-(x_vals + sqrt(a/b)).^2/(2*a/(2*b)));
pdf_unstable = (1/sqrt(2*pi)) * exp(-x_vals.^2/(2/(a)));
total_pdf = pdf_stable1/(2*sqrt(a/(2*b))) + pdf_stable2/(2*sqrt(a/(2*b))) + pdf_unstable/(2*sqrt(1/a));plot(x_vals, total_pdf, 'r', 'LineWidth', 2);
title('状态分布与理论概率密度');
xlabel('状态变量 x');
ylabel('概率密度');
legend('仿真分布', '理论分布');
grid on;%% 8. 频谱分析
% 计算FFT
Fs = 1/dt; % 采样频率
NFFT = 2^nextpow2(N); % FFT点数
f = Fs/2*linspace(0,1,NFFT/2+1); % 频率轴% 输入信号频谱
signal_fft = fft(signal_norm, NFFT)/N;
signal_psd = 2*abs(signal_fft(1:NFFT/2+1));% 输出信号频谱
output_fft = fft(output, NFFT)/N;
output_psd = 2*abs(output_fft(1:NFFT/2+1));% 噪声频谱
noise_fft = fft(noise_norm, NFFT)/N;
noise_psd = 2*abs(noise_fft(1:NFFT/2+1));% 绘制频谱
figure('Name', '频谱分析', 'Position', [100, 100, 1200, 800]);
subplot(3,1,1);
plot(f, 10*log10(signal_psd), 'b', 'LineWidth', 1.5);
title('输入信号频谱');
xlabel('频率 (Hz)');
ylabel('功率谱密度 (dB)');
grid on;
xlim([0, 0.2]);subplot(3,1,2);
plot(f, 10*log10(noise_psd), 'r', 'LineWidth', 1.5);
title('噪声频谱');
xlabel('频率 (Hz)');
ylabel('功率谱密度 (dB)');
grid on;
xlim([0, 0.2]);subplot(3,1,3);
plot(f, 10*log10(output_psd), 'k', 'LineWidth', 1.5);
hold on;
plot([f0, f0], ylim, 'r--', 'LineWidth', 1.5); % 标记信号频率
title('输出信号频谱');
xlabel('频率 (Hz)');
ylabel('功率谱密度 (dB)');
legend('输出频谱', '信号频率');
grid on;
xlim([0, 0.2]);%% 9. 信噪比分析
% 计算信噪比改善
signal_power = sum(signal_norm.^2);
noise_power_input = sum(noise_norm.^2);
input_snr = 10*log10(signal_power / noise_power_input);output_signal_power = sum(signal_norm.^2); % 假设信号成分不变
output_noise_power = sum((output - signal_norm).^2);
output_snr = 10*log10(output_signal_power / output_noise_power);snr_improvement = output_snr - input_snr;fprintf('输入信噪比: %.2f dB\n', input_snr);
fprintf('输出信噪比: %.2f dB\n', output_snr);
fprintf('信噪比改善: %.2f dB\n', snr_improvement);%% 10. 随机共振参数扫描
% 扫描噪声强度
D_values = linspace(0.01, 0.5, 20);
snr_values = zeros(size(D_values));for i = 1:length(D_values)% 生成新噪声temp_D = D_values(i);switch noise_typecase 'gaussian'temp_noise = sqrt(2*temp_D*dt) * randn(1, N);case 'uniform'temp_noise = sqrt(12*temp_D*dt) * (rand(1, N) - 0.5);case 'colored'raw_noise = sqrt(2*temp_D*dt) * randn(1, N);alpha = 0.9;temp_colored_noise = zeros(1, N);temp_colored_noise(1) = raw_noise(1);for j = 2:Ntemp_colored_noise(j) = alpha*temp_colored_noise(j-1) + (1-alpha)*raw_noise(j);endtemp_noise = temp_colored_noise;end% 归一化噪声if normalize_outputswitch normalization_methodcase 'minmax'temp_noise_norm = (temp_noise - min(temp_noise)) / (max(temp_noise) - min(temp_noise));case 'zscore'temp_noise_norm = (temp_noise - mean(temp_noise)) / std(temp_noise);endelsetemp_noise_norm = temp_noise;end% 运行系统temp_x = zeros(1, N);temp_x(1) = 0.1;for j = 1:N-1dx = (a*temp_x(j) - b*temp_x(j)^3 + signal_norm(j) + temp_noise_norm(j)) * dt;temp_x(j+1) = temp_x(j) + dx;end% 计算输出信噪比temp_output = temp_x;if normalize_outputswitch normalization_methodcase 'minmax'temp_output = (temp_x - min(temp_x)) / (max(temp_x) - min(temp_x));case 'zscore'temp_output = (temp_x - mean(temp_x)) / std(temp_x);endendtemp_output_signal_power = sum(signal_norm.^2);temp_output_noise_power = sum((temp_output - signal_norm).^2);if temp_output_noise_power > 0snr_values(i) = 10*log10(temp_output_signal_power / temp_output_noise_power);elsesnr_values(i) = inf;end
end% 绘制信噪比曲线
figure('Name', '随机共振曲线', 'Position', [100, 100, 800, 600]);
plot(D_values, snr_values, 'bo-', 'LineWidth', 1.5, 'MarkerSize', 8);
xlabel('噪声强度 D');
ylabel('输出信噪比 (dB)');
title('随机共振曲线');
grid on;% 标记最佳噪声强度
[~, idx_max] = max(snr_values);
hold on;
plot(D_values(idx_max), snr_values(idx_max), 'rs', 'MarkerSize', 12, 'MarkerFaceColor', 'r');
text(D_values(idx_max), snr_values(idx_max), ...sprintf(' 最佳D=%.3f\n SNR=%.2fdB', D_values(idx_max), snr_values(idx_max)), ...'FontSize', 10, 'VerticalAlignment', 'bottom');%% 11. 参数敏感性分析
% 扫描信号频率
f0_values = linspace(0.01, 0.1, 10);
snr_freq = zeros(size(f0_values));for i = 1:length(f0_values)temp_f0 = f0_values(i);temp_signal = A * sin(2*pi*temp_f0*t + phase);if normalize_outputswitch normalization_methodcase 'minmax'temp_signal_norm = (temp_signal - min(temp_signal)) / (max(temp_signal) - min(temp_signal));case 'zscore'temp_signal_norm = (temp_signal - mean(temp_signal)) / std(temp_signal);endelsetemp_signal_norm = temp_signal;end% 运行系统temp_x = zeros(1, N);temp_x(1) = 0.1;for j = 1:N-1dx = (a*temp_x(j) - b*temp_x(j)^3 + temp_signal_norm(j) + noise_norm(j)) * dt;temp_x(j+1) = temp_x(j) + dx;end% 计算输出信噪比temp_output = temp_x;if normalize_outputswitch normalization_methodcase 'minmax'temp_output = (temp_x - min(temp_x)) / (max(temp_x) - min(temp_x));case 'zscore'temp_output = (temp_x - mean(temp_x)) / std(temp_x);endendtemp_output_signal_power = sum(temp_signal_norm.^2);temp_output_noise_power = sum((temp_output - temp_signal_norm).^2);if temp_output_noise_power > 0snr_freq(i) = 10*log10(temp_output_signal_power / temp_output_noise_power);elsesnr_freq(i) = inf;end
end% 绘制频率响应曲线
figure('Name', '频率响应', 'Position', [100, 100, 800, 600]);
plot(f0_values, snr_freq, 'go-', 'LineWidth', 1.5, 'MarkerSize', 8);
xlabel('信号频率 f_0 (Hz)');
ylabel('输出信噪比 (dB)');
title('频率响应特性');
grid on;%% 12. 归一化方法比较
% 比较不同归一化方法
methods = {'none', 'minmax', 'zscore'};
results = zeros(length(methods), 3); % [SNR, 执行时间, 状态方差]for m = 1:length(methods)method = methods{m};tic; % 开始计时% 设置归一化选项if strcmp(method, 'none')normalize_output = false;elsenormalize_output = true;normalization_method = method;end% 重新运行仿真% ...（此处省略重复代码，实际应用中应封装为函数）% 简化：使用之前的结果if strcmp(method, 'none')snr_val = output_snr;var_val = var(output);elseif strcmp(method, 'minmax')snr_val = output_snr; % 假设相同var_val = var(output); % 假设相同elsesnr_val = output_snr; % 假设相同var_val = var(output); % 假设相同endexec_time = toc; % 结束计时results(m, :) = [snr_val, exec_time, var_val];
end% 显示比较结果
fprintf('\n归一化方法比较:\n');
fprintf('方法\t\tSNR(dB)\t执行时间(s)\t状态方差\n');
for m = 1:length(methods)fprintf('%s\t%.2f\t%.4f\t\t%.4f\n', methods{m}, results(m,1), results(m,2), results(m,3));
end%% 13. 高级分析：Kramers逃逸速率
% 计算理论逃逸速率
V1 = -a/(2*b); % 势阱深度
kramers_rate = (1/(2*pi)) * sqrt(a^2/(2*b)) * exp(-2*V1/(D));fprintf('\n理论Kramers逃逸速率: %.4f Hz\n', kramers_rate);% 从仿真数据估计逃逸速率
% 通过状态穿越中点次数估计
crossings = 0;
threshold = 0; % 中点位置
for i = 2:length(output)if output(i-1) < threshold && output(i) > thresholdcrossings = crossings + 1;elseif output(i-1) > threshold && output(i) < thresholdcrossings = crossings + 1;end
endestimated_rate = crossings / (2 * t_end); % 双向穿越，除以2
fprintf('估计逃逸速率: %.4f Hz\n', estimated_rate);%% 14. 保存结果
save('stochastic_resonance_results.mat', 't', 'signal', 'noise', 'output', 'D_values', 'snr_values');
随机共振原理与实现
1. 基本概念
随机共振是一种非线性现象，当弱周期性信号与适当强度的噪声共同作用时，系统输出信噪比会显著提高。这种现象最早在古气象学中被发现，现已广泛应用于物理学、生物学和神经科学等领域。
2. 双稳态系统模型
本程序使用一维双势阱系统作为随机共振的基本模型：
V(x) = -a*x²/2 + b*x⁴/4
系统动力学方程：
dx/dt = a*x - b*x³ + A*sin(2πf₀t) + Γ(t)
其中：
x是状态变量A*sin(2πf₀t)是输入的弱周期信号Γ(t)是高斯白噪声
3. 归一化处理
程序实现了两种归一化方法：
-
最小-最大归一化：将信号映射到[0,1]区间
x_norm = (x - min(x)) / (max(x) - min(x)) -
Z-score标准化：使信号均值为0，标准差为1
x_norm = (x - mean(x)) / std(x)
4. 数值求解方法
使用Euler-Maruyama方法求解随机微分方程：
x(t+Δt) = x(t) + [a*x(t) - b*x(t)³ + s(t)]Δt + √(2DΔt)·ξ(t)
其中ξ(t)是标准正态分布随机数。
关键功能模块
1. 信号生成模块
% 周期性信号
signal = A * sin(2*pi*f0*t + phase);% 高斯白噪声
noise = sqrt(2*D*dt) * randn(1, N);
2. 归一化处理模块
switch normalization_methodcase 'minmax'signal_norm = (signal - min(signal)) / (max(signal) - min(signal));case 'zscore'signal_norm = (signal - mean(signal)) / std(signal);
end
3. 系统动力学模块
for i = 1:N-1dx = (a*x(i) - b*x(i)^3 + signal_norm(i) + noise_norm(i)) * dt;x(i+1) = x(i) + dx;
end
4. 信噪比分析模块
signal_power = sum(signal_norm.^2);
noise_power = sum((output - signal_norm).^2);
snr = 10*log10(signal_power / noise_power);
5. 参数扫描模块
for i = 1:length(D_values)% 更新噪声强度temp_noise = sqrt(2*D_values(i)*dt) * randn(1, N);% 运行系统% ...% 计算SNR% ...
end
仿真结果分析
1. 时域信号分析
- 输入信号：显示弱周期信号
- 噪声：显示添加的随机噪声
- 系统输出：显示经过双稳态系统处理后的信号
2. 相空间分析
- 绘制状态变量x与其导数的关系图
- 展示系统在相空间中的运动轨迹
3. 状态分布分析
- 比较仿真得到的状态分布与理论概率密度函数
- 验证双稳态系统的两个势阱和势垒
4. 频谱分析
- 显示输入信号、噪声和系统输出的频谱
- 标记信号频率位置，观察共振峰
5. 随机共振曲线
- 扫描噪声强度D，计算对应的输出信噪比
- 确定最佳噪声强度（随机共振点）
6. 参数敏感性分析
- 扫描信号频率f₀，分析系统频率响应特性
- 确定系统的最佳工作频率范围
参考代码 归一化随机共振代码 www.youwenfan.com/contentcnm/83391.html
扩展功能与应用
1. 多稳态系统扩展
% 三稳态系统
V = @(x) -a*x.^2/2 + b*x.^4/4 - c*x.^6/6;
dVdx = @(x) -a*x + b*x.^3 - c*x.^5;
2. 耦合系统实现
% 两个耦合的双稳态系统
dx1 = (a*x1 - b*x1^3 + coupling*(x2-x1) + signal + noise1)*dt;
dx2 = (a*x2 - b*x2^3 + coupling*(x1-x2) + signal + noise2)*dt;
3. 自适应随机共振
% 自适应调整噪声强度
if output_snr < target_snrD = D * 1.1; % 增加噪声
elseD = D * 0.9; % 减少噪声
end
4. 生物医学信号处理
% 处理ECG信号
ecg_signal = load('ecg_data.mat');
A = 0.2; % 弱信号幅值
D = 0.15; % 最优噪声强度
5. 图像处理应用
% 图像去噪
image_noisy = imread('noisy_image.png');
for each pixel% 应用随机共振去噪
end
实际应用场景
1. 微弱信号检测
- 地震波检测
- 心电信号分析
- 机械故障诊断
2. 神经科学
- 神经元放电分析
- 听觉感知研究
- 视觉信息处理
3. 量子系统
- 量子比特读出
- 约瑟夫森结电路
- 激光系统
4. 金融数据分析
- 股票价格预测
- 汇率波动分析
- 市场周期检测
结论
本MATLAB实现提供了完整的归一化随机共振仿真框架：
- 精确建模：实现双稳态系统的随机微分方程
- 灵活归一化：支持多种归一化方法
- 全面分析：包含时域、频域、相空间分析
- 参数扫描：支持噪声强度和信号频率扫描
- 性能评估：计算信噪比改善和逃逸速率
通过调整参数和扩展模型，该程序可用于：
- 微弱信号检测系统优化
- 非线性系统特性研究
- 生物医学信号处理
- 金融数据分析
