Two-variable boosting bifurcation in a hyperchaotic map and its hardware implementation
对这篇论文做笔记和复现以及分析
1.Abstract:
1.对离散图中混沌序列的振荡幅度进行无损调整的报道较少，本文设计了一个新的基于三角函数的二维超混沌的模型，
2.通过改变初始变量和系统参数展示了偏移增强分岔和偏移增强共存吸引子，而且还显示了两个状态变量相对于二维图中任意参数的偏移增强
3.引入新的控制参数也可以使超混沌映射表现出丰富的多稳态现象
4.再通过硬件测试和PRNG来验证
2.二维超混沌映射的数学模型
xn+1 = xn + a sin(yn)
yn+1 = b sin(xn)
求出其不动点和其特征值
特征方程：P(λ) = (λ − 1)λ − ab cos(μπ )
P1(λ) = (λ − 1)λ − ab
P2(λ) = (λ − 1)λ + ab
由u是否为奇偶数来决定
再分析其SE(谱熵复杂度) 谱熵通常用于描述信号的复杂性和不规则性。
给出计算SE的函数
function SE=SEShannon(x)%不需要修改这个程序
%函数名称：SEShannon
%函数功能：计算序列的谱熵复杂度
%输入参数x；x为混沌序列
%输出参数： SE；SE为输出的混沌序列复杂度
N=length(x);
flag=0;
x=x-mean(x);
for i=1:N
if x(i)~=0
flag=1;
end
end
if flag==0
SE=0;
return;
end
Y=fft(x);
% if Y(1)==0
% Y(1)=1e-30;
% end
Xk=(abs(Y).^2)./N;
Xk=Xk(1:(floor(N/2)));
ptot=sum(Xk);
PK=Xk./ptot;
P=0;
for i=1:N/2
if PK(i)~=0
P=P-PK(i)*log(PK(i));
end
end
se2=sum(P);
SE=se2/log(N/2);
此为这个代码的具体解释
这个MATLAB函数用于计算输入混沌序列的谱熵复杂度。谱熵通常用于描述信号的复杂性和不规则性。下面是代码的详细解释：
1. `function SE=SEShannon(x)`：这一行定义了MATLAB函数，函数名为`SEShannon`，输入参数是混沌序列`x`，输出参数是复杂度值`SE`。
2. `N=length(x);`：获取输入序列`x`的长度，这将用于后续计算中。
3. `flag=0;`：定义一个名为`flag`的标志，初始化为0，稍后用于检查序列`x`是否全为零。
4. `x=x-mean(x);`：将序列`x`减去其均值，以确保序列的均值为零。
5. `for i=1:N`：开始一个循环，遍历序列`x`中的每个元素。
6. `if x(i)~=0`：检查序列中的元素是否不等于零。
7. `flag=1;`：如果有任何一个元素不等于零，将`flag`标志设置为1。
8. `end`：结束条件判断。
9. `if flag==0`：检查`flag`标志是否仍然为0，如果是，则表示序列`x`中的所有元素都为零，因此谱熵复杂度为0。
10. `SE=0;`：在这种情况下，将`SE`设置为0，表示复杂度。
11. `return;`：返回结果并退出函数。
12. `end`：结束`if flag==0`条件判断。
13. `Y=fft(x);`：对序列`x`进行傅里叶变换，得到频域表示。
14. `Xk=(abs(Y).^2)./N;`：计算傅里叶变换后的频谱功率密度，即每个频率分量的平方幅值除以`N`。
15. `Xk=Xk(1:(floor(N/2)));`：截取频谱的前一半，因为对称性质，后一半是多余的。
16. `ptot=sum(Xk);`：计算所有频率分量的总功率。
17. `PK=Xk./ptot;`：计算每个频率分量的功率占比。
18. `P=0;`：初始化一个变量`P`，用于存储计算谱熵的中间值。
19. `for i=1:N/2`：开始一个循环，遍历频率分量。
20. `if PK(i)~=0`：检查功率占比是否不等于零，以避免log(0)。
21. `P=P-PK(i)*log(PK(i));`：根据Shannon熵的计算公式，计算每个频率分量的熵贡献，并将其累加到`P`中。
22. `end`：结束条件判断。
23. `se2=sum(P);`：将所有频率分量的熵贡献相加，得到`se2`。
24. `SE=se2/log(N/2);`：计算谱熵复杂度，将`se2`除以log(N/2)，其中`N`是输入序列的长度的一半。这个值即为谱熵的复杂度。
这个函数主要用于计算给定混沌序列的谱熵，谱熵用于描述序列的复杂性和不规则性。
此为主函数：
clear;clc;
L=100;
C=linspace(-3,3,L);%横坐标区间
C1=linspace(-3,3,L);%纵坐标区间
SE=zeros(L,L);
dot=50000;
for k1=1:L
b=C1(k1);%纵坐标
for k=1:L
a=C(k); %横坐标
x(1)=0.1;
y(1)=0.1;
for n=1:dot
x(n+1)=x(n)+a*sin(y(n));
y(n+1)=b*sin(x(n));
end
SE(k1, k)=SEShannon(x(10000:end));
clear x y
end
disp(k1)
end
figure
colormap(flipud(hot));
contourf(C,C1,SE, 'LineStyle','none');
xlabel('\ita')
ylabel('\itb')
set(gca,'FontName','Times New Roman','box','on','linewidth',1,'fontsize',16);
toc
3.基于初始值的初始提升分叉和提升共存吸引子
首先是令初值(x0, y0) = (0.1, 0.1)和系统的参数a=1.8
代码如下：
为分岔图
clc
clear all
Q=eye(2);
dot=100000;
b=-3:0.01:3; %b
for k1=1:length(b)
m=b(k1);%分叉参数
x(1)=0.1;
y(1)=0.1;
a=1.8;
for n=1:dot
x(n+1)=x(n)+a*sin(y(n));
y(n+1)=m*sin(x(n));
end
B1(k1,:)=x(end-100:end);
clear x y z
end
figure(8)
plot(b,B1,'r.','Markersize',4)
xlabel('\ita','fontsize',16)
ylabel('\itx_n','fontsize',16)
set(gca,'FontName','Times New Roman','linewidth',1,'fontsize',16);
hold on
用y为纵坐标修改B1(k1,:)=x(end-100:end);即可
李指数图：
clc
clear all
Q=eye(2);
dot=40000;
b=2;
b=-8:0.01:8; %分岔参数范围
% % u=-0.5;d=1.3;m=-0.6;
% u=0.8;r=0.05;d=-0.56;k=0.99;
% d=-0.64;c=1;
for k1=1:length(b)
k=b(k1);%参数
x(1)=0.1;
y(1)=0.1;
a=1.8;
% z(1)=0.1;%初值
for n=1:dot
% x(n+1)=x(n)+a*x(n)*sin(y(n))+b;
% y(n+1)=y(n)+x(n); %超混沌映射系统
% x(n+1)=1+y(n)-a*x(n)^2;
% y(n+1)=b*x(n);henon映射
%
% x(n+1)=x(n)-k*x(n)*cos(y(n));
% y(n+1)=y(n)+x(n);
x(n+1)=x(n)+a*sin(y(n));
y(n+1)=k*sin(x(n));%ND文章的模型
% x(n+1)=k*d*sin(c*pi*x(n))+(((2*y(n)-1)^2)*x(n));
% y(n+1)=x(n)+y(n);%Chaos文章的模型
% x(n+1)=k*(d*(2*y(n)-1)^2*x(n)+u*sin(pi*x(n))+m*cos(z(n))*x(n));
% y(n+1)=y(n)+x(n);%基于忆阻的混沌映射
% z(n+1)=x(n)+z(n);
% x(n+1)=x(n)*u*((y(n)^2-1)+z(n))+d*sin(pi*x(n)); no.4
% y(n+1)=y(n)+k*x(n);
% z(n+1)=c*cos(r*n+r);NO.4论文设想
% x(n+1)=k*d*sin(c*pi*x(n))+(((2*y(n)-1)^2)*x(n));
% y(n+1)=x(n)+y(n);%Chaos文章的模型
Xx=1;
Xy=a*cos(y(n));
Yx=k*cos(x(n));
Yy=0;
J=[Xx Xy ;Yx Yy ;];
B=J*Q;
[Q,R]=qr(B);
L1(n+1)=log(abs(R(1,1)));
L2(n+1)=log(abs(R(2,2)));
end
B1(k1,:)=x(end-100:end);
L11(k1)=sum(L1(1000:end))/30000;
L22(k1)=sum(L2(1000:end))/30000;
end
% figure%分岔图
% plot(A,B1,'b.','Markersize',4)
% % xlim([0.2 1.25])
% % ylim([-1.5 2.5])
% xlabel('\itb','fontsize',16)
% ylabel('\itx\rm(\itn\rm)','fontsize',16)
% set(gca,'FontName','Times New Roman','linewidth',1,'fontsize',16);
figure%李指图
plot(b,L11,'r','linewidth',2);
grid on
hold on
plot(b,L22,'b','linewidth',2);
% grid on
% hold on
% plot(A,L33,'k','linewidth',2);
% xlim([0.2 1.25])
% ylim([-5 1.1])
xlabel('\itb','fontsize',16)
ylabel('LEs','fontsize',16)
legend('\lambda_1','\lambda_2','fontsize',16)
set(gca,'FontName','Times New Roman','linewidth',1,'fontsize',16);
接下来的用a为行坐标的分岔图和李指数同上进行适当修改即可
对于改变参数和初始值的对称，绘制吸引子和迭代序列
tic;
a=1.8;b=2;%参数
dot=50000;
x(1)=0.1;%初值
y(1)=0.1；
for n=1:dot
x(n+1)=x(n)+asin(y(n));
y(n+1)=bsin(x(n));
end
figure(22)
plot(x(1:400),’’)
xlabel(‘\itn’,’fontsize’,16)
ylabel(‘\itx\rm(\itn\rm)’,’fontsize’,16)
set(gca,’FontName’,’Times New Roman’,’linewidth’,1,’fontsize’,16);
hold on
x(1)=-0.1;%初值
y(1)=-0.1;
for n=1:dot
x(n+1)=x(n)+a*sin(y(n));
y(n+1)=b*sin(x(n));
end
plot(x(1:400),'')
plot(x(1000:15000),y(1000:15000),'r.')
xlabel('\itx','fontsize',18)
ylabel('\ity','fontsize',18)
set(gca,'FontName','Times New Roman','linewidth',1,'fontsize',16);
figure(49)
plot(x(1000:15000),y(1000:15000),'b.')
xlabel('\itx','fontsize',18)
ylabel('\ity','fontsize',18)
set(gca,'FontName','Times New Roman','linewidth',1,'fontsize',16);
hold on
只列举了其中a=1.8,b=2的例子，其他参数类似
上述结果表示：本文提出的二维超混沌映射具有复杂的动态特性，其动态行为与两个系统参数密切相关
0—1测试的介绍：
公式(10)中，M(n) 表示平均平方位移，表示随时间 n 变化的均方根位移。它的计算基于维度函数 p(n) 和 q(n)。
维度函数 p(n) 和 q(n) 是在复平面上的均值分量，它们的计算涉及傅里叶变换，详见前面的讨论。
Mc(n)（公式(11)）是平均平方位移 M(n) 的计算，它是根据时间序列 x(n) 的差的平方来计算的。它是在序列上计算的，其中 j 表示不同的时间偏移，N 表示总的采样点数。
E(φ)（公式(12)）是时间序列 x(n) 的平均值。
Kc（公式(13)）是 M(n) 的渐进增长率。它表示 M(n) 随时间 n 变化的趋势。Kc 的计算是取 M(n) 的对数，然后除以 log(n)，其中 n 表示时间。
这些公式的推导需要一些数学上的分析和极限过程，通常通过极限分析来获得。总体来说，这些公式用于描述混沌序列的平方位移和随时间的增长趋势，以便判断序列的混沌性质。如果 Kc 接近 0 或 1，那么序列更可能是混沌的。
0-1检测方法的思路是将混沌序列映射到复平面上，并根据分布情况来判断序列的混沌性。下面是更详细的步骤：
获取混沌序列：首先，你有一个混沌序列，假设这个序列是x(n)，n表示序列的离散时间。
分别计算p(n)和q(n)：为了将序列映射到复平面，你需要计算序列在不同方向上的均值分量。
p(n) 表示序列的实部在不同方向上的均值分量。
q(n) 表示序列的虚部在不同方向上的均值分量。
这些分量是通过傅里叶变换的实部和虚部来计算的。
映射到复平面：使用p(n)作为x坐标，q(n)作为y坐标，你将得到一系列点在复平面上的分布。
分析复平面上的点的分布：观察复平面上的点的分布情况，特别是它们是否靠近单位圆（与混沌性相关）。通过计算逐渐增长率Kc（参考步骤8的公式），你可以根据这些点的分布来判断混沌性。
判断混沌性：根据步骤4中复平面上的点的分布情况和计算的Kc值，你可以判断混沌序列的特性。具体的判断标准通常是Kc接近0和1，越接近0或1表示序列越可能是混沌的。
结果：根据Kc的值和复平面上的点的分布，你可以判断输入的序列是否具有混沌性质。
总之，0-1检测方法通过将混沌序列映射到复平面上，分析复平面上的点的分布来判断混沌性。p(n)和q(n)分别表示序列的实部和虚部在不同方向上的均值分量，是检测的关键步骤之一
这里代码较长就不列举代码了（已打包）
这个测试方法也为检测混沌行为的一种方法（记忆）
下面介绍 基于初始值的初始提升分叉
通过改变初始值(x0, y0)=(0.1+2k π ,0.1)(k=±2, ±1, 0)
and (x0, y0)=(−0.1 + 2k π , −0.1)(k=±2, ±1, 0)
其正负值的分岔图是对称的
而且按+2kπ的方式改变初始值 得到的一些指标相同，这表明混沌映射具有良好的初始可控性和鲁棒性。
其分岔图和吸引子的代码与前文的类似，就不再列举，只需改变初值即可
介绍一下PE(排列熵)
对于计算排列熵（Permutation Entropy，PE），它是一种用来衡量时间序列复杂性的方法，特别是在混沌理论和动力系统研究中应用广泛。它通过考虑时间序列中数据点的排列来量化时间序列的规则性或可预测性。以下是如何计算排列熵的步骤：
选择嵌入维度：首先，您需要确定一个适当的嵌入维度，用“m”表示。这个维度与您用于构建排列的数据点数相关。
创建长度为“m”的向量：将您的时间序列数据进行处理，通过在时间序列上滑动一个窗口来创建一组长度为“m”的向量。窗口每次移动一个时间步长，每个向量包含连续的“m”个数据点。
计算排列：对于每个长度为“m”的向量，计算数据点的排列方式数，以量化数据点的排列顺序。
计算概率：计算每种唯一排列在您的向量中出现的概率。概率是通过将特定排列出现的次数除以排列的总数来确定的。
计算排列熵：最后，您可以计算排列熵，通常表示为PE(m)，通过将所有唯一排列的概率相加并取负自然对数。PE的通常公式如下：
PE(m) = -∑[P(π) * log(P(π))]
其中：
PE(m) 表示嵌入维度“m”的排列熵。
P(π) 表示特定排列“π”出现的概率。
您可以使用不同的嵌入维度进行此计算，比如PE(3)和PE(5)，它们分别表示嵌入维度为3和5时的排列熵。排列熵可以测量时间序列的复杂性，PE值较低表示更规律和可预测的行为，而PE值较高表示更大的复杂性和不可预测性。
代码：
function [PE3, PE5] = calculatePermutationEntropy(data)
% 输入参数 data：包含时间序列数据的向量
% 设置嵌入维度
m3 = 3; % 嵌入维度为3
m5 = 5; % 嵌入维度为5
% 初始化排列熵
PE3 = 0;
PE5 = 0;
% 数据长度
N = length(data);
% 计算排列熵 PE(3)
for i = 1:(N - m3 + 1)
subsequence = data(i:(i + m3 - 1));
[~, order] = sort(subsequence);
PE3 = PE3 - log(factorial(m3)) + log(factorial(m3) - sum(diff(order) == 1));
end
PE3 = PE3 / (N - m3 + 1);
% 计算排列熵 PE(5)
for i = 1:(N - m5 + 1)
subsequence = data(i:(i + m5 - 1));
[~, order] = sort(subsequence);
PE5 = PE5 - log(factorial(m5)) + log(factorial(m5) - sum(diff(order) == 1));
end
PE5 = PE5 / (N - m5 + 1);
end
4.1/4.2由参数控制的增长的分岔图和吸引子和迭代序列
xn+1 = xn + a sin(yn)
yn+1 = b sin(xn) + mπ
新加入参数m
分岔图与吸引子：
周期性的改变m，分岔图依然呈现出boost分岔图
通过m的周期性的改变 控制参数可以无损地控制混沌/超混沌迭代序列的动态幅度，
其代码只要改变其数学模型即可
4.3 当混沌映射的参数设定为a = 1.93，b = 2.83时，控制参数m的周期性切换可以产生有趣的吸引子共存现象
说明控制参数可以有效地控制吸引子迭代序列的动态幅值
当控制参数不是整数时 显示了吸引子的非对称共存
4.4 2D超混沌映射序列的性能分析
将本文提出的混沌映射和现有混沌映射的动态性能( 2015 ~ 2022年)进行对比，
用LLE,SE的方法比较出本文的混沌映射有更好的动力学性能
5.1硬件的实现
基于stm32和ad数模转换器来实现
将二维混沌映射烧录到单片机中，用示波器双通道显示其混沌序列
由于手上没有开发板就不演示了
表明了本文超混沌映射硬件实现的正确性和可行性，并表明硬件设备中混沌映射的初始控制混沌序列和三参数混沌序列具有较好的鲁棒性。
这里解释一下鲁棒性：
在上下文中，”鲁棒性”（robustness）是指系统或现象对外部干扰、变化或扰动的抵抗能力和稳定性。在文中提到 “initial controllability of the hyperchaotic map is robust”，意味着超混沌地图（hyperchaotic map）在不同的初始条件或干扰情况下能够保持其可控性和稳定性。
这意味着即使超混沌地图经历一些不同的初始条件、参数变化或其他扰动，它仍然能够被有效地控制和保持其预期的行为
这种鲁棒性对于许多应用是非常重要的，因为它确保了系统的可控性和可靠性，即使在不确定性或外部变化的情况下也是如此
5.2在伪随机数发生器中的应用
将混沌序列转换为二进制序列通常涉及到浮点数到二进制的转换。在文中提到使用IEEE 754浮点标准，下面是将浮点数转换为二进制的一般步骤：
获取浮点数：首先，你有一个混沌序列，例如X(n) = {x0, x1, x2, ...}。
规范化：确保你的浮点数在0和1之间，因为IEEE 754标准通常用于规范化浮点数。如果你的混沌序列范围不在0到1之间，你可能需要将其缩放到这个范围内。
IEEE 754标准：该标准定义了如何将浮点数表示为二进制数。它通常包括三个部分：符号位、指数位和尾数位。符号位用于表示正数或负数，指数位用于表示浮点数的阶码，尾数位包含实际的数值。
转换：按照IEEE 754标准的要求，将规范化的浮点数转换为二进制。这通常涉及到将小数部分转换为二进制，并调整指数部分以保持浮点数的值。
填充：确保你的二进制序列长度足够，以满足你的需要。文中提到二进制流的长度不少于106，因此，你需要将二进制序列填充到至少10^6位
代码
function binary = floatToBinary(number)
% 检查输入是否为非负数
if number < 0
error('输入必须是非负数.');
end
% IEEE 754标准的位数
bitCount = 32; % 单精度浮点数，你可以根据需要进行调整
% 初始化二进制字符串
binary = '';
% 处理特殊情况
if number == 0
binary = '0';
return;
end
% 处理符号位
binary = '0'; % 假设输入为正数
% 将浮点数转换为二进制表示
while number > 0 && length(binary) < bitCount
number = number * 2;
if number >= 1
binary = [binary '1'];
number = number - 1;
else
binary = [binary '0'];
end
end
% 如果二进制字符串太短，进行填充
while length(binary) < bitCount
binary = [binary '0'];
end
end
再使用NIST套件来检测伪随机数生成器（PRNG）输出的随机性和质量
按照使用说明输入序列，即可得到的相应P-T值
这表明所提出的超混沌映射能够产生具有高随机性的伪随机数。
个人总结：
本文提出了一种简单的2D混沌映射模型，通过改变初始值实现boost分岔
再通过引入新的参数，也可已实现状态变量的boost分岔
而且可以通过初始状态和控制参数的周期性变化来调整实现状态变量的无损调整