对一个轴承故障信号Acos(2Πft)采用随机共振进行降噪，对非线性方程用四阶龙格库塔输出，其参数a=1 b=1 噪声D=0.5 写出matlab代码

以下是一种可能的 Matlab 代码实现：

% 定义参数
f = 100; % 轴承故障信号频率
a = 1; % 非线性方程参数a
b = 1; % 非线性方程参数b
D = 0.5; % 噪声强度

% 生成轴承故障信号
t = 0:0.001:1; % 时间向量
x = cos(2*pi*f*t); % 原始信号

% 加入高斯白噪声
y = x + D*randn(size(x));

% 预处理随机共振参数
Q = 10; % 质量因子
f0 = f; % 共振频率等于信号频率
w0 = 2*pi*f0; % 角频率
K = w0^2/Q; % 力学增益
C = 2*Q*w0; % 阻尼系数

% 定义非线性方程
f_nl = @(t, x) [x(2); -K*x(1) - C*x(2) - a*x(1)^3 - b*x(1)*x(2)^2];

% 利用四阶龙格库塔求解非线性方程
[t, z] = ode45(f_nl, t, [0;0]);

% 显示结果
figure;
subplot(211);
plot(t, x, 'b', t, y, 'g');
xlabel('Time (s)');
ylabel('Signal');
title('Original Signal (blue) and Noisy Signal (green)');
legend('Original', 'Noisy');
subplot(212);
plot(t, z(:,1), 'r');
xlabel('Time (s)');
ylabel('Signal');
title('Filtered Signal');
legend('Filtered');

% 随机共振滤波函数
function y = sr_filter(x, Q, f0, D)
w0 = 2*pi*f0; % 角频率
K = w0^2/Q; % 力学增益
C = 2*Q*w0; % 阻尼系数
dt = 0.001; % 采样周期
t = 0:dt:(length(x)-1)*dt; % 时间向量
xi = randn(size(x)); % 高斯白噪声
y = zeros(size(x)); % 输出向量
for n = 3:length(x)
y(n) = -K*y(n-1) - C*y(n-2) + x(n) - D*xi(n);
end
end

% 生成高斯白噪声函数
function y = randn(size)
y = sqrt(-2*log(rand(size))) .* cos(2*pi*rand(size));
end