随机共振的Matlab源代码，供大家一起讨论。

2012-10-19 10:58| 发布者: aspen| 查看: 1319| 评论: 0|原作者: wuliping1|来自: 振动论坛

摘要: 下面是我修改别人的随机共振系统的Matlab代码，供大家讨论，本人也是初步学习，有关随机共振的详细内容，希望和大家一起学习，共同进步^_^ 主程序： clear all clc fs=5; %参样频率 f=0.01; %信号频率 Ts=1/ ...

下面是我修改别人的随机共振系统的Matlab代码，供大家讨论，本人也是初步学习，有关随机共振的详细内容，希望和大家一起学习，共同进步^_^
主程序：

clear all
clc
fs=5; %参样频率
f=0.01; %信号频率
Ts=1/fs;%参样时间
h=1/fs; %时间步长
t=0:Ts:4095*Ts;
D=0.6; %噪声强度
%双稳态系统参数
a=1;
b=1;
s=0.3*sin(2*pi*0.01*t); %信号
x1=s+sqrt(2*D)*randn(size(t)); %噪声
%输入无噪信号傅立叶变换
y=fft(s,4096);
pyy=y.*conj(y)/4096;
ff=fs*(0:2048)/4096;
figure(1);
subplot(2,1,1);plot(t,s);
title('输入无噪信号');xlabel('时间t/s');ylim([-0.5,0.5]);ylabel('信号幅度A');
subplot(2,1,2);plot(ff,pyy(1:2049));
xlabel('频率f/Hz');ylabel('频谱幅度');xlim([0,0.05]);title('输入无噪信号的频谱');
%输入加噪信号傅立叶变换
y=fft(x1,4096);
pyy=y.*conj(y)/4096;
ff=fs*(0:2048)/4096;
figure(2)
subplot(2,1,1);plot(t,x1);
title('输入加噪噪信号');xlabel('时间t/s');ylabel('信号幅度A');
subplot(2,1,2);plot(ff,pyy(1:2049));
xlabel('频率f/Hz');ylabel('频谱幅度');xlim([0,0.05]);ylim([0,1500]);title('输入加噪信号的频谱');
%四阶龙格库塔法对双稳态输出信号求解
x=sr(a,b,h,x1);
%输出信号求傅立叶变换
y=fft(x,4096);
py=y.*conj(y)/4096;
ff=fs*(0:2048)/4096;
figure(3);
subplot(2,1,1);plot(t,x);
title('输出信号');xlabel('时间t/s');ylabel('信号幅度A');
subplot(2,1,2);plot(ff,py(1:2049));
xlabel('频率f/Hz');ylabel('频谱幅度');xlim([0,0.05]);ylim([0,1500]);title('输出信号的频谱');

子程序：解四阶龙格库塔法
function x=sr(a,b,h,x1)
x=zeros(1,length(x1));
for i=1:length(x1)-1
k1=h*(a*x(i)-b*x(i).^3+x1(i));
k2=h*(a*(x(i)+k1/2)-b*(x(i)+k1/2).^3+x1(i));
k3=h*(a*(x(i)+k2/2)-b*(x(i)+k2/2).^3+x1(i+1));
k4=h*(a*(x(i)+k3)-b*(x(i)+k3).^3+x1(i+1));
x(i+1)=x(i)+(1/6)*(k1+2*k2+2*k3+k4);
end

收藏分享邀请

账号		自动登录	找回密码
密码			我要加入

随机共振的Matlab源代码，供大家一起讨论。

最新评论

相关分类

下级分类