xcsyb 发表于 2006-6-22 09:39

[求助]求随机共振源程序

谢谢各位!!!

风花雪月 发表于 2006-6-25 09:19

回复:(xcsyb)[求助]求随机共振源程序

随机共振的基本含义是指一个非线性双稳或多稳系统,当在小周期信号驱动下不足使系统的输出在两个稳态之间跳跃,只能在一个稳态上。而在弱噪声和小周期调制信号共同作用下,随着输入噪声强度的增加,系统输出功率谱中调制信号的频率处出现一个峰值.当噪声强度达到某一适当值时,系统在两定态之间进行跃迁,输出信号的峰值达到最大——“共振”。以后随着噪声强度的增强,输出信号的峰值逐渐减小。总结产生随机共振现象的条件,大致可以归纳为以下三个要素: <BR>(1)非线性系统要存在两个(或多个)稳态。 <BR>(2)输入信号,并且输入信号还不足引起系统在不同稳态之间跃迁。 <BR>(3)噪声。它既可以是系统外部的随机驱动力,也可以来自系统内部相互作用引起的涨落。 <BR>

风花雪月 发表于 2006-6-25 09:19

回复:(xcsyb)[求助]求随机共振源程序

下面是matlab写的一个随机共振的例子,你可以参考一下%programshijin.m
x=-2:0.1:2;区间范围
b=1;势阱参数,决定势阱的深度与宽度
c=1;

figure(1)
subplot(2,1,1);分区作图
y=-b.*x.^2/2+c.*x.^4/4;势函数
plot(x,y);作二维图形
subplot(2,1,2);
b=3:0.1:7;给定b,c的变化区间
c=2:0.1:6;
h=moviein(40);影片动画
fori=1:40;
y=-b(i).*x.^2/2+c(i).*x.^4/4;
plot(x,y);
h(:,i)=getframe;
end
title('势函数随b,c的变化情况')坐标标注
subplot(2,1,1)
title('无信号及噪声输入时系统的势函数')
xlabel('x');
ylabel('U(x)');
subplot(2,1,2)
xlabel('x');
ylabel('U(x)');
movie(h,1);
figure(2)
b=1;给定b,c加入周期调制信号
c=1;
w=pi;
z=0;
subplot(2,1,1)信号较小时
A=2;
m=moviein(20);影片动画
=ode23('xinhao1',,,[],b,z,c,A,w);ode命令解微分方程
n=length(t);
fori=1:fix(n/2)
axis([-22-510]);
y=-b.*x.^2/2+c.*x.^4/4-A.*x.*sin(w*t(i));
plot(x,y);
holdon
uu=yy(i,2)*yy(i,2)/2-b.*yy(i,1).^2/2+c.*yy(i,1).^4/4-A.*yy(i,1).*sin(w*t(i));
uu=-b.*yy(i,1).^2/2+c.*yy(i,1).^4/4-A.*yy(i,1).*sin(w*t(i));
plot(yy(i,1),uu,'.r','markersize',10)
holdoff
m(:,i)=getframe;
end
title('加入周期信号较小,不足以翻转')
xlabel('x');
ylabel('U(x)');
v=moviein(20);
subplot(2,1,2)周期信号较大时
A=3;
=ode23('xinhao1',,,[],b,z,c,A,w);
n=length(t);
fori=1:fix(n/2)
y=-b.*x.^2/2+c.*x.^4/4-A.*x.*sin(w*t(i));
plot(x,y);
holdon
uu=yy(i,2)*yy(i,2)/2-b.*yy(i,1).^2/2+c.*yy(i,1).^4/4-A.*yy(i,1).*sin(w*t(i));
uu=-b.*yy(i,1).^2/2+c.*yy(i,1).^4/4-A.*yy(i,1).*sin(w*t(i));
plot(yy(i,1),uu,'.r','markersize',10)
holdoff
v(:,i)=getframe;
end
title('加入周期调制信号较大,可以翻转')
holdon
xlabel('x');
ylabel('U(x)');
axis([-22-510])
holdon
movie(v)
t=0:0.1:10;
n=length(t);
figure(3)加入信号同时加入随即噪声
A=2;
b=1;
m=moviein(20);
fori=1:fix(n/2)
r=2*rand(1)随即量
y=-b.*x.^2/2+c.*x.^4/4-A.*x.*sin(w*t(i))-r.*x;
plot(x,y);
holdon
=oula('f',0,5,rand(1),0.1,r);调用欧拉函数解方程
u=-b.*s(i).^2/2+c.*s(i).^4/4-A.*s.*sin(w*t(i))-r*s(i);
plot(s(i),u(i),'.r','markersize',10);
holdoff
m(:,i)=getframe;
end
title('加入随机噪声的情况,参数选取同figure(2)中未翻转时');
xlabel('x');
ylabel('U(x)');
movie(m,1)
figure(4)波形图及频谱分析
x=-2:0.1:2;
w=50*pi;
b=1;
c=1
A=1;
t=0:0.01:10;
n=length(t);
subplot(2,1,1);
fori=1:n
y(i)=b.*x(1)-c.*x(1).^3+A.*cos(w*t(i))+0.5*rand(1);
end
plot(t,y)
title('信号的波形图');
xlabel('t');
ylabel('U(x)');
subplot(2,1,2)
z=fft(y);
z(1)=[];
n=length(z);
p=10^(-5).*abs(z(1:n/2)).^2
m=1/2;
f=100*(1:n/2)/(n/2)*m;
plot(f,p)
title('信号的频谱图,w=50*pi,f=25');


%oula.m子程序:欧拉法解方程
function=oula(fun,x0,xf,y0,h,r)
n=fix((xf-x0)/h);
x(1)=x0;
y(1)=y0;
x(n)=0;
y(n)=0;
fori=1:(n-1)
x(i+1)=x0+i*h;
y1=y(i)+h*feval(fun,x(i),y(i),r);调用方程
y2=y(i)+h*feval(fun,x(i+1),y1,r);
y(i+1)=(y1+y2)/2;
end


%f.m方程文件
%参数b,c,A,w为1,故省略。
functionf=f(t,x,r)
f=x-x.^3+cos(t)-r;


%odefilexinhao1.mode文件
functionydot=xinhao1(t,x,flag,b,z,c,A,w)
ydot=[b*x(1)-c*x(1)^3+A*cos(w*t)+z;...
b*x(2)-3*c*x(1)*x(1)*x(2)-A*w*sin(w*t)];

[ 本帖最后由 风花雪月 于 2008-4-30 10:30 编辑 ]

xcsyb 发表于 2006-6-29 22:41

回复:(风花雪月)回复:(xcsyb)[求助]求随机共振源...

谢谢!!!

nanlaibeiqu 发表于 2008-10-15 17:32

回复 板凳 风花雪月 的帖子

为什么我运行不了楼主的程序呢,提示
??? Error: File: sr1.m Line: 123 Column: 1
Function definitions are not permitted at the prompt or in scripts.

风花雪月 发表于 2008-10-28 17:18

原帖由 nanlaibeiqu 于 2008-10-15 17:32 发表 http://www.chinavib.com/forum/images/common/back.gif
为什么我运行不了楼主的程序呢,提示
??? Error: File: sr1.m Line: 123 Column: 1
Function definitions are not permitted at the prompt or in scripts.

都不知道你看过我上面的帖子没有,自己仔细看看吧,而不是盲目的粘贴运行

lihaitao123 发表于 2013-4-3 21:43

风花雪月 发表于 2008-10-28 17:18 static/image/common/back.gif
都不知道你看过我上面的帖子没有,自己仔细看看吧,而不是盲目的粘贴运行

谢谢风花雪月,我刚接触随机振动,太有帮助了。。

lordkibun 发表于 2013-6-5 22:42

风花雪月 发表于 2006-6-25 09:19 static/image/common/back.gif
下面是matlab写的一个随机共振的例子,你可以参考一下

[ 本帖最后由 风花雪月 于 2008-4-30 10:30 编辑 ]

您好。我也是在看这段程序。想问下,您用的是哪个版本的matlab?
页: [1]
查看完整版本: [求助]求随机共振源程序