采用Hilbert变换求耦合系统相位差的问题

liuhongjuan

如题，我采用如下程序，不管C值取什么，得到的结果都是一条直线，而论文中给出的结果如附件所示，请问大侠们，这是怎么回事，多谢！
程序：
function y=rossler3(t,x)
global omiga
global C
y=[-(1+omiga)*x(2)-x(3)+C*(x(4)-x(1));
    (1+omiga)*x(1)+0.15*x(2);
     0.2+x(3)*(x(1)-10);
   -(1-omiga)*x(5)-x(6)+C*(x(1)-x(4));
    (1-omiga)*x(4)+0.15*x(5);
     0.2+x(6)*(x(4)-10)];

clear all
global omiga
global C
omiga=0.015;
C=input('input the value of controling parameter c =');
% 微分方程的初始值。
x0=[0.15,1.23,1.09,2.0,1.25,1.76];
tspan=[0,2000];
options=odeset('RelTol',1e-10);
[t,y]=ode45(@rossler3,tspan,x0,options);
% 希尔伯特变换相位。
hy=hilbert(y);
%[fnor,t]=instfreq(hy);
%plot(t,fnor*fs);
hy=angle(hy);
hy=unwrap(hy);
% x1
phi1=atan(imag(hy(:,1))./y(:,1));
% x2
phi2=atan(imag(hy(:,4))./y(:,4));
% 相位差。
phi=phi1 - phi2;
plot(t,phi,'r')