%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1 混沌系统
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
混沌系统具有对初始条件或者系统参数敏感的特性。当系统参数选取不同的值时，
系统能够表现出倍周期分叉(Bifurcation diagram) 和混沌 (Chaos) 等复杂现象。
下面主要介绍一下用三阶微分方程表示的连续混沌系统。
（1）Lorenz 系统
dx/dt=a*(y-x)
dy/dt=b*x-y-x*z
dz/dt=x*y-c*z
其中a,b,c是系统参数。当a=10, b=28, c=8/3 时，Lorenz 系统是混沌的。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
tic;h=0.01;
x=zeros(1,10001);x(1)=1;
y=zeros(1,10001);y(1)=1;
z=zeros(1,10001);z(1)=1;
for r=1:10000
x(r+1)=x(r)+h*(10*(y(r)-x(r)));
y(r+1)=y(r)+h*(-x(r)*z(r)+28*x(r)-y(r));
z(r+1)=z(r)+h*(x(r)*y(r)-8/3*z(r));
end
t=toc;disp([Time for run is:,num2str(t)]);
plot(x,y,k.,marker,.,markersize,1)
title(洛伦兹方程)
print(-dbmp,洛伦兹方程时间序列)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
(2) Chen 系统
dx/dt=a*(y-x)
dy/dt=(c-a)*x+c*y-x*z
dz/dt=x*y-b*z
其中a,b,c是系统参数。当a=35, b=3, c=28 时系统是混沌的。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
tic;h=0.0001;a=35;b=3;c=28;
x=zeros(1,100001);%x(1)=1;
y=zeros(1,100001);y(1)=1;
z=zeros(1,100001);%z(1)=1;
for r=1:100000
x(r+1)=x(r)+h*a*(y(r)-x(r));
y(r+1)=y(r)+h*((c-a)*x(r)+c*y(r)-x(r)*z(r));
z(r+1)=z(r)+h*(x(r)*y(r)-b*z(r));
end
t=toc;disp([Time for run is:,num2str(t)]);
plot(x(10000:100001),y(10000:100001),k.,marker,.,markersize,1)
title(Chen 系统)
print(-dbmp,Chen 系统)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
(3)Chua 系统
dx/dt=-apha*x+apha*y-apha*f
dy/dt=x-y+z
dz/dt=beta*z
其中f=b*x+(a-b)*(|y+1|-|y-1|)/2, apha,beta,a,b是系统参数。
当参数apha=10, beta=14.87, a=-1.27, b=-0.68时，Chua系统是混沌的。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
tic;h=0.0001;apha=10;beta=14.87;a=-1.27;b=-0.68;
x=zeros(1,10001);x(1)=1;
y=zeros(1,10001);y(1)=1;
z=zeros(1,10001);z(1)=1;
for r=1:10000
f(r)=b*x(r)+(a-b)*(abs(y(r)+1)-abs(y(r)-1))/2;
x(r+1)=x(r)+h*(-apha*x(r)+apha*y(r)-apha*f(r));
y(r+1)=y(r)+h*(x(r)-y(r)+z(r));
z(r+1)=z(r)+h*(beta*z(r));
end
t=toc;disp([Time for run is:,num2str(t)]);
plot(x,y,k.,marker,.,markersize,1)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
(4) Rossler 系统
dx/dt=-y-z
dy/dt=x+a*y
dz/dt=c-b*z+x*z
其中a,b,c是系统参数。当a=0.2, b=5.7, c=0.2 时， Rossler 系统是混沌的。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
tic;h=0.01;
x=zeros(1,10001);x(1)=1;
y=zeros(1,10001);y(1)=1;
z=zeros(1,10001);z(1)=1;
for r=1:10000
x(r+1)=x(r)+h*(-y(r)-z(r));
y(r+1)=y(r)+h*(x(r)+0.2*y(r));
z(r+1)=z(r)+h*(0.2-5.7*z(r)+x(r)*z(r));
end
t=toc;disp([Time for run is:,num2str(t)]);
plot(x,y,k.,marker,.,markersize,1)
title(Rossler系统)
print(-dbmp, Rossler系统xz)
%%%%%%%%%%%%%%%%%四阶Rossler振荡器
clear
tic;h=0.01;
x=zeros(1,10001);x(1)=1;
y=zeros(1,10001);y(1)=1;
z=zeros(1,10001);z(1)=1;
w=zeros(1,10001);w(1)=1;
for r=1:10000
x(r+1)=x(r)+h*(-y(r)-z(r));
y(r+1)=y(r)+h*(x(r)+0.25*y(r)+w(r));
z(r+1)=z(r)+h*(3+x(r)*z(r));
w(r+1)=z(r)+h*(-0.5*z(r)+0.05*w(r));
end
t=toc;disp([Time for run is:,num2str(t)]);
plot(x,y,k.,marker,.,markersize,1)
title(Rossler振荡器)
print(-dbmp,Rossler振荡器xy)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
(5) van der Pol 振子
dx/dt=x-x^3/3-y+p+q*cos(w*t)
dy/dt=c*(x+a-b*y)
其中a,b,c,p,q,w是系统参数。当a=0.7, b=0.8, c=0.1, p=0.0, q=0.74, w=1.0 时，van der Pol 是混沌的。
除此之外，著名的连续混沌还有Duffing 振子和 Lu (吕) 系统。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
tic;h=0.01;a=0.7;b=0.8;c=0.1;p=0.0;q=0.74;w=1.0;
x=zeros(1,10001);x(1)=1;
y=zeros(1,10001);y(1)=1;
z=zeros(1,10001);z(1)=1;
for r=1:10000
x(r+1)=x(r)+h*(x(r)-x(r)^3/3-y(r)+p+q*cos(w*h));
y(r+1)=y(r)+h*c*(x(r)+a-b*y(r));
end
t=toc;disp([Time for run is:,num2str(t)]);
plot(x,y,k.,marker,.,markersize,1)
title(van der Pol 振子)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
(6) Lu (吕) 系统
clear
tic;h=0.001;a=36;b=3;c=20;
x=zeros(1,100001);x(1)=1;
y=zeros(1,100001);y(1)=1;
z=zeros(1,100001);z(1)=1;
for r=1:100000
x(r+1)=x(r)+h*a*(y(r)-x(r));
y(r+1)=y(r)+h*(-x(r)*z(r)+c*y(r));
z(r+1)=z(r)+h*(x(r)*y(r)-b*z(r));
end
t=toc;disp([Time for run is:,num2str(t)]);
plot3(x(10000:100001),y(10000:100001),z(10000:100001),k.,marker,.,markersize,1)
title(Lu (吕) 系统)
print(-dbmp, Lu (吕) 系统)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2 gray变换
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%gray变换
clear
tic
T=[1,0,0,0,0,0,0,0;
1,1,0,0,0,0,0,0;
0,1,1,0,0,0,0,0;
0,0,1,1,0,0,0,0;
0,0,0,1,1,0,0,0;
0,0,0,0,1,1,0,0;
0,0,0,0,0,1,1,0;
0,0,0,0,0,0,1,1];
G=zeros(1,256);
r=0:255;
g=de2bi(r,left-msb,8);
G(r+1)=[128,64,32,16,8,4,2,1]*rem(T*g(r+1,:),2);
toc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
help de2bi %可以利用de2bi进行改进
用2进制gray进行256*256图象加密
clear
A=imread(lena.bmp);
T=[1,0,0,0,0,0,0,0;
1,1,0,0,0,0,0,0;
0,1,1,0,0,0,0,0;
0,0,1,1,0,0,0,0;
0,0,0,1,1,0,0,0;
0,0,0,0,1,1,0,0;
0,0,0,0,0,1,1,0;
0,0,0,0,0,0,1,1];
for i=1:256;
for j=1:256;
t=double(A(i,j));
u7=fix(t/128);
u6=fix((t-128*u7)/64);
u5=fix((t-128*u7-64*u6)/32);
u4=fix((t-128*u7-64*u6-32*u5)/16);
u3=fix((t-128*u7-64*u6-32*u5-16*u4)/8);
u2=fix((t-128*u7-64*u6-32*u5-16*u4-8*u3)/4);
u1=fix((t-128*u7-64*u6-32*u5-16*u4-8*u3-4*u2)/2);
u0=t-128*u7-64*u6-32*u5-16*u4-8*u3-4*u2-2*u1;
u=[u7,u6,u5,u4,u3,u2,u1,u0];
B(i,j)=[128,64,32,16,8,4,2,1]*rem(T^2*u,2);
end
end
imshow(uint8(B))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 利用128*128幻方对图像进行加密 （骑士游加密方法与之类似）
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%按顺序找出矩阵中元素的坐标：
clear
n=128;
a=magic(n);
[i,j,v]=find(a);%v为从第一列开始a中的元素，i,j为其位置坐标
i=i;
j=j;
v=v;
[E,index]=sort(v);%index为E的索引。即：index(t)个元素为E(t)
for t=1:n^2
y1(t)=i(index(t));
y2(t)=j(index(t));
end
A=imread(camera.gif);
[m,n]=size(A);
B=zeros(m,n);
for N=1:2
for t=1:m^2-1;
B(y1(t+1),y2(t+1))=A(y1(t),y2(t));
end
B(y1(1),y2(1))=A(y1(m^2),y2(n^2));
A=B;
end
imshow(uint8(B));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4 用面包师变换对256*256图象进行1次加密
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
A=imread(lena.bmp);
[M,N]=size(A);
for i=1:M
for j=1:N
if (1<=i)&(i<=M/2)&(1<=j)&(j<=N/2)
x=2*i;
y=2*j-1;
elseif (1<=i)&(i<=M/2)&(N/2 x=2*i-1;
y=2*N-2*j-2;
elseif (M/2 x=2*i-M;
y=2*j;
else
x=2*i-M-1;
y=2*N-2*j+1;
end
B(257-y,x)=A(257-j,i);
end
end
imshow(B);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5 对图像进行SVD变换%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
[A,map]=imread(lena.bmp); %读取图像
%A=X(120:220,120:220); %将图像的一部分取出，并命名为A
[nx,ny]=size(A); %算出A矩阵的大小??
I=eye(nx,nx); %取一个和A同样大小的单位矩阵
for i=nx:-1:nx/4 %将单位矩阵I的部分主对角线元素设为0
I(i,i)=0; %以便和A 的特征值矩阵相乘后???
end %消去一部分的特征值?
[u d v]=svd(double(A)); %对A矩阵进行SVD变换
B=u*I*d*v;
B=uint8(B); %将 A 矩阵还原回来并设为 B
subplot(221);imshow(A,map) %将原始图像显示出来
subplot(222);imhist(A)
subplot(223);imshow(B,map) %将近似图像显示出来
subplot(224);imhist(B)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
6 用Arnold变换对256*256灰度图象进行n次加密
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Arnoldtran(image,N)
A=imread(image);
[m,n]=size(A);
%tic;
for i=0:m-1
for j=0:n-1
Y=Arnoldtransf(i,j,N,m);
B(m-Y(1),Y(2)+1)=A(m-i,j+1);
end
end
%t=toc;
%disp([Time for run is:,num2str(t)]);
subplot(121);imshow(A);title(原始图像 );
subplot(122);imshow(B);title(加密图像);
%imwrite(B,192.bmp);
function Y=Arnoldtransf(x,y,N,m)
Y=[x,y];
for k=1:N
t=x;
x=mod(x+y,m);
y=mod(t+2*y,m);
end
Y=[x,y];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
周期
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#include
void main()
{
cout<<"求二维Arnold变换在数字图象矩阵的维数从2到点512时的周期"< for(int i=2;i<=1024;i++)
{
int x=i-1;
int y=i-1;
int j=0;
int t;
do
{
t=x;
x=(7*x+10*y)%i;
y=(5*t+7*y)%i;
j=j+1;
}
while((x!=(i-1))||(y!=(i-1)));
cout< }
}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
7 此程序为Z形扫描的逆过程，由一位好朋友提供，在此表示感谢！
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function a=izigzag(v,m,n)
a=zeros(m,n);
a(1,1)=v(1);
counter=2;
for k=3:1:m+n
b1=bound1(k,a);
b2=bound2(k,a);
if mod(k,2)==1
for p=b1:1:b2
a(p,k-p)=v(counter);
counter=counter+1;
end
else
for p=b2:-1:b1
a(p,k-p)=v(counter);
counter=counter+1;
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% inner function
function b1=bound1(k,a)
[m,n]=size(a);
for b1=1:1:m
q=k-b1;
if q<=n
return
end
end
%%%%%%%%%%%%%%%%%%%%%%
function b2=bound2(k,a)
[m,n]=size(a);
for q=1:1:n
b2=k-q;
if b2<=m
return
end
end
%%%%%%%%%%%%%%%%%%%%%
