T=cputime;
h=0.02;
t=0:h:100-h;
fs=1/h;
N=length(t);
x(1)=0;
A=0.3;
f=0.02;
Y=[];
q=1;
P1=[];
for d=0.1:0.1:1
for n=2:N xA=0;yA=0;xB=0;yB=0;
for j=1:n-1
if j==1,a=(n-1)^(q+1)-(n-q-1)*n^q;
else a=(n-j+2)^(q+1)+(n-j)^(q+1)-2*((n-j+1)^(q+1));
end
b=((h^q)/q)*((n-j+1)^q-(n-j)^q);
u(j)=x(j)-x(j)^3+A*cos(2*pi*f*h*j)+sqrt(d)*randn(1);
C(j)=a*u(j);
D(j)=b*u(j);
%xA=sum(C1);
%yA=sum(C2);
xA=xA+C(j);
xB=xB+D(j);
end
xp(n)=x(1)+1/gamma(q)*xB;
x(n)=x(1)+h^q/gamma(2+q)*(xp(n)-xp(n)^3+A*cos(2*pi*f*h*n)+sqrt(d)*randn(1))+h^q/gamma(2+q)*xA;
end
y0=xcorr(x,'unbiased');
N=length(y0);
y0=y0/y0((N-1)/2+1);
%归一化
y=fft(y0);
y1=abs(y);
N1=round(f*(N-1)/fs+1);
SNout=y1(N1)/((sum(y1)-y1(N1))/length(y1));
P1=[P1,SNout];
end
d=0.1:0.1:1;
plot(d,P1,'o');
cputime=T;
一直都是BUSY，一直繁忙，不知道什么时候能运行出来，是不是代码有问题啊？