MATLAB归一化随机共振代码

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. 生成噪声switchnoise_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);fori=2:Ncolored_noise(i)=alpha*colored_noise(i-1)+(1-alpha)*raw_noise(i);endnoise=colored_noise;otherwiseerror('未知噪声类型');end%% 4. 归一化处理% 归一化输入信号ifnormalize_outputswitchnormalization_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('未知归一化方法');endelsesignal_norm=signal;noise_norm=noise;end%% 5. 双稳态系统模型% 初始化状态变量x=zeros(1,N);x(1)=0.1;% 初始条件% 系统动力学方程% dx/dt = a*x - b*x^3 + signal + noisefori=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. 归一化输出ifnormalize_outputswitchnormalization_methodcase'minmax'output=(x-min(x))/(max(x)-min(x));case'zscore'output=(x-mean(x))/std(x);otherwiseoutput=x;endelseoutput=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. 频谱分析% 计算FFTFs=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));fori=1:length(D_values)% 生成新噪声temp_D=D_values(i);switchnoise_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);forj=2:Ntemp_colored_noise(j)=alpha*temp_colored_noise(j-1)+(1-alpha)*raw_noise(j);endtemp_noise=temp_colored_noise;end% 归一化噪声ifnormalize_outputswitchnormalization_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;forj=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;ifnormalize_outputswitchnormalization_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);iftemp_output_noise_power>0snr_values(i)=10*log10(temp_output_signal_power/temp_output_noise_power);elsesnr_values(i)=inf;endend% 绘制信噪比曲线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));fori=1:length(f0_values)temp_f0=f0_values(i);temp_signal=A*sin(2*pi*temp_f0*t+phase);ifnormalize_outputswitchnormalization_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;forj=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;ifnormalize_outputswitchnormalization_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);iftemp_output_noise_power>0snr_freq(i)=10*log10(temp_output_signal_power/temp_output_noise_power);elsesnr_freq(i)=inf;endend% 绘制频率响应曲线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, 执行时间, 状态方差]form=1:length(methods)method=methods{m};tic;% 开始计时% 设置归一化选项ifstrcmp(method,'none')normalize_output=false;elsenormalize_output=true;normalization_method=method;end% 重新运行仿真% ...（此处省略重复代码，实际应用中应封装为函数）% 简化：使用之前的结果ifstrcmp(method,'none')snr_val=output_snr;var_val=var(output);elseifstrcmp(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');form=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;% 中点位置fori=2:length(output)ifoutput(i-1)<threshold&&output(i)>threshold crossings=crossings+1;elseifoutput(i-1)>threshold&&output(i)<threshold crossings=crossings+1;endendestimated_rate=crossings/(2*t_end);% 双向穿越，除以2fprintf('估计逃逸速率: %.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. 归一化处理

程序实现了两种归一化方法：

1. 最小-最大归一化：将信号映射到[0,1]区间

   x_norm = (x - min(x)) / (max(x) - min(x))

2. 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. 归一化处理模块

switchnormalization_methodcase'minmax'signal_norm=(signal-min(signal))/(max(signal)-min(signal));case'zscore'signal_norm=(signal-mean(signal))/std(signal);end

3. 系统动力学模块

fori=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. 参数扫描模块

fori=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/contentcsn/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. 自适应随机共振

% 自适应调整噪声强度ifoutput_snr<target_snr D=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');foreach pixel% 应用随机共振去噪end

实际应用场景

1. 微弱信号检测

• 地震波检测
• 心电信号分析
• 机械故障诊断

2. 神经科学

• 神经元放电分析
• 听觉感知研究
• 视觉信息处理

3. 量子系统

• 量子比特读出
• 约瑟夫森结电路
• 激光系统

4. 金融数据分析

• 股票价格预测
• 汇率波动分析
• 市场周期检测

结论

本MATLAB实现提供了完整的归一化随机共振仿真框架：

1. 精确建模：实现双稳态系统的随机微分方程
2. 灵活归一化：支持多种归一化方法
3. 全面分析：包含时域、频域、相空间分析
4. 参数扫描：支持噪声强度和信号频率扫描
5. 性能评估：计算信噪比改善和逃逸速率

通过调整参数和扩展模型，该程序可用于：

• 微弱信号检测系统优化
• 非线性系统特性研究
• 生物医学信号处理
• 金融数据分析