用FPGA设计超混沌数字电路

WC_scanf

已于 2025-03-03 10:23:13 修改

阅读量434

点赞数 6

分类专栏： FPGA设计文章标签： fpga开发

于 2025-02-27 13:42:32 首次发布

本文链接：https://blog.youkuaiyun.com/A_scanf_/article/details/142783022

版权

FPGA设计专栏收录该内容

17 篇文章

订阅专栏

介绍：

超混沌系统的复杂动态特性使其在信号加密中表现出色，能有效防止破解。同时与模拟电路相比，数字电路对噪声不敏感，稳定性更高，随着技术进步，超混沌数字电路在保密通信、信息安全及复杂系统控制等领域展现出巨大潜力，成为当前科研热点之一。

当前，实现混沌系统主要有两种方法：

第一种是离散算法编写底层硬件代码如，欧拉算法、龙格库塔等，实现更加灵活。

第二种则是借助 Matlab 的 DSP BUILDER 库搭建混沌数字电路，无需深入编程，但使用并不灵活，容易受限于硬件资源的限制。

本文介绍运用改良欧拉法实现超混沌系统的过程，并成功实现了一个具备隐藏多稳态吸引子的混沌系统发生器。

原理：

超混沌Lorenz系统是对经典Lorenz系统的扩展，它是由爱德华·诺伦·洛伦兹（Edward Norton Lorenz）发展而来的一种非线性动力学系统。爱德华·洛伦兹是一位美国气象学家，于1963年提出了经典Lorenz系统，用于描述大气运动中的热对流现象。经典Lorenz系统是一种简化的三维模型，由三个非线性常微分方程组成。然而，洛伦兹后来意识到，通过对经典Lorenz系统进行一些修改，可以得到更加复杂和混沌的动力学行为。于是，他引入了一些非线性项和额外的驱动力，从而得到了超混沌Lorenz系统。其数学公式如下

$\left\{\begin{array}{l}\dot{x}=a(y-x)+w \\ \dot{y}=c x-y-x z \\ \dot{z}=x y-b z \\ w=-y z+r w\end{array}\right.$

选取改良欧拉法对超混沌Lorenz系统进行离散化。其数学公式如下

$\begin{cases} \tilde{y}_{n + 1} = y_n + hf(t_n, y_n) \\ y_{n + 1} = y_n + \frac{h}{2}\left[f(t_n, y_n)+f(t_{n + 1}, \tilde{y}_{n + 1})\right] \end{cases}$

实验目的：用FPGA实现数字化的超混沌Lorenz电路，并通过设置初始值来影响混沌状态

（参考链接：基于FPGA的混沌系统实现_fpga混沌序列是怎样产生的-优快云博客）

结果

模块框图：

模块代码

hyperchaos模块代码

module hyperchaos(
	input wire clk,
	input wire rst_n,
	input wire [15:0] ran_num,

	output  signed [31:0] x,
	output  signed [31:0] y,
	output  signed [31:0] z,
	output  signed [31:0] w
);

    parameter t = 12; //影响算法计算值
	
	parameter S0=3'b000,
		      S1=3'b001,
		      S2=3'b010,
		      S3=3'b100;

	parameter a  = 8   ,  //2.4.8.16 /1 2 3 4 
		      b  = 2   ,  //2.4 /1 2
		      c   = 32'b0_00010_10101010101010101010101011;
		
	reg signed [31:0] x_nt;
	reg signed [31:0] y_nt;
	reg signed [31:0] z_nt;
	reg signed [31:0] w_nt;
 

    reg signed [31:0] x_n_temp   ;//计算的中间值，用来存储x_n的值
    reg signed [31:0] y_n_temp   ;
    reg signed [31:0] z_n_temp   ;
    reg signed [31:0] w_n_temp   ;
			  
    reg signed [31:0] dx;
	reg signed [31:0] dy;
	reg signed [31:0] dz;
	reg signed [31:0] dw;

    reg signed [31:0] dx_temp;
	reg signed [31:0] dy_temp;
	reg signed [31:0] dz_temp;
	reg signed [31:0] dw_temp;

	wire signed [63:0]xz_64;//两个32位的数相乘是64位
	wire signed [63:0]xy_64;
	wire signed [63:0]cz_64;
	wire signed [63:0]yz_64;

	assign xz_64 = x_nt * z_nt;//乘法
	assign xy_64 = x_nt * y_nt;
	assign yz_64 = y_nt * z_nt;
	assign cz_64 = c * z_nt;

	wire signed [31:0]xz;//对xz_64位数进行截取后的结果
	wire signed [31:0]xy;
	wire signed [31:0]cz;
	wire signed [31:0]yz;

    wire signed [31:0]bx;

    assign bx =  (x_nt <<< 5) - (x_nt <<< 2); //乘28


	assign xz = {xz_64[63],xz_64[56:26]} ;//截取规则为保留符号位，摒弃低26位小数位（可忽略）和高6位整数位（一般做了压缩后，都是为0）
	assign xy = {xy_64[63],xy_64[56:26]} ;
	assign yz = {yz_64[63],yz_64[56:26]} ;
	assign cz = {cz_64[63],cz_64[56:26]} ;

    reg signed [31:0] xz_extend_6;//比例压缩，由于怕计算过程内部出现比五位整数位大的数，导致内部计算溢出，x,y,z都压缩64。例如dx = xy +z,将x,y,z压缩64，原=原式变换为d64x = 64x*64y +64z等价于dx=32xy+z;即n项式要乘于32^（n-1） 
    reg signed [31:0] xy_extend_6;//同上
    reg signed [31:0] yz_extend_6;

	reg [2:0]state;
    reg flag;
    reg flag_p;
    
	
	always@(posedge clk or negedge rst_n) begin
		if(!rst_n) begin
			state <= S0;
			end
		else case (state) 
			S0 : state <= S1; 
			S1 : begin
                if(flag == 1) begin
                    state <= S3;
                end
                else begin
                    state <= S2; 
                end
            end
			S2 : state <= S0;
			S3 : state <= S0;
			default: state <= S0;			
		endcase
    end 

	//设置传入值
	always@(posedge clk or negedge rst_n) begin
		if(!rst_n) begin
			x_nt <= {16'b0_00000_0000010000,ran_num};
			y_nt <= {16'b0_00000_0000001000,ran_num};
			z_nt <= {16'b0_00000_0000001000,ran_num};
			w_nt <= {16'b0_00000_0000001000,ran_num};
		end
		else if(state == S2) begin //根据k1 更新传入值 来计算新k2
            x_nt <= x_nt + (dx>>>(t-1));//对x_n进行关于算法计算
            y_nt <= y_nt + (dy>>>(t-1));
            z_nt <= z_nt + (dz>>>(t-1));   
            w_nt <= w_nt + (dw>>>(t-1));   
		end
		else if(state == S3) begin //根据k1 k2更新传入值 来计算新的微分值
			x_nt <= x_n_temp + (dx>>>t) + (dx_temp>>>t);
			y_nt <= y_n_temp + (dy>>>t) + (dy_temp>>>t);
			z_nt <= z_n_temp + (dz>>>t) + (dz_temp>>>t);
			w_nt <= w_n_temp + (dw>>>t) + (dw_temp>>>t);
		end
	end		

    //flag 用于记录是否已经计算
    always@(posedge clk or negedge rst_n) begin
		if(!rst_n) begin
            flag <= 0; 
		end
        else if(state == S1) begin
            flag <= 0;
		end
		else if(state == S2) begin
            flag <= 1;
		end
	end

	

	always@(posedge clk or negedge rst_n) begin
		if(!rst_n) begin
            flag_p <= 1; 
		end
        else if(state == S0) begin
            flag_p <= 0;
		end
		else if(state == S3) begin
            flag_p <= 1;
		end
	end

    //保存x_n的值
    always@(posedge clk or negedge rst_n) begin
		if(!rst_n) begin
            x_n_temp <= 0;
            y_n_temp <= 0;
            z_n_temp <= 0;
            w_n_temp <= 0;
		end
		else if(state == S2) begin
            x_n_temp <= x_nt;
            y_n_temp <= y_nt;
            z_n_temp <= z_nt;
            w_n_temp <= w_nt;
		end
	end

    //保存dx的值
    always@(posedge clk or negedge rst_n) begin
		if(!rst_n) begin
            dx_temp <= 0;
            dy_temp <= 0;
            dz_temp <= 0;
            dw_temp <= 0;
		end
		else if(state == S2) begin
            dx_temp <= dx;
            dy_temp <= dy;
            dz_temp <= dz;
            dw_temp <= dw;
		end
	end

    //乘64
    always@(posedge clk or negedge rst_n) begin
		if(!rst_n) begin
            xz_extend_6 <= 0;
            xy_extend_6 <= 0;
            yz_extend_6 <= 0;
		end
		else if(state == S0) begin
			xz_extend_6 <= xz <<< 6;
            xy_extend_6 <= xy <<< 6;
            yz_extend_6 <= yz <<< 6;
		end
	end


	//计算参数
	always@(posedge clk or negedge rst_n) begin
		if(!rst_n) begin
			dx <= 0;
		end
		else if(state == S1) begin
			dx <= (y_nt -x_nt) * a + w_nt;
		end
	end

	always@(posedge clk or negedge rst_n) begin
		if(!rst_n) begin
			dy <= 0;
		end
		else if(state == S1) begin
			dy <= bx - xz_extend_6 - y_nt;
		end
	end

	always@(posedge clk or negedge rst_n) begin
		if(!rst_n) begin
			dz <= 0;
		end
		else if(state == S1) begin
			dz <= xy_extend_6 -cz;
		end
	end

	always@(posedge clk or negedge rst_n) begin
		if(!rst_n) begin
			dw <= 0;
		end
		else if(state == S1) begin
			dw <= - yz_extend_6 - (((w_nt << 1) + w_nt) >> 1); 
		end
	end

    wire signed [31:0]x_temp;//结果，由于可能为负数，一般加上一个正数后才输出，保证输出全为正数
	wire signed [31:0]y_temp;
	wire signed [31:0]z_temp;
    wire signed [31:0]w_temp;

    assign x_temp= (flag_p)?x_nt:x_temp;
    assign y_temp= (flag_p)?y_nt:y_temp;
    assign z_temp= (flag_p)?z_nt:z_temp;	
    assign w_temp= (flag_p)?w_nt:w_temp;
   
    assign x = x_temp+ 32'b0000_0100_0000_0000_0000_0000_0000_0000;//加个1保证输出为正
    assign y = y_temp+ 32'b0000_0100_0000_0000_0000_0000_0000_0000;
    assign z = z_temp+ 32'b0000_0100_0000_0000_0000_0000_0000_0000;
    assign w = w_temp+ 32'b0000_0100_0000_0000_0000_0000_0000_0000;

endmodule

仿真验证部分

hyperchaos_tb模块代码

//tb文件
`timescale 1ns/1ns
`define clock_period 20

module hyperchaos_tb;

	reg clk;
	reg rst_n;
	wire [31:0]x,y,z,w;

hyperchaos hyperchaos_inst  
(
	.clk    (clk    ),
	.rst_n  (rst_n  ),
    .ran_num(16'h5555),

	.x      (x      ),
	.y      (y      ),
	.z      (z      ),
    .w      (w      )
);

integer handle_x,handle_y,handle_z,handle_w;
integer i;
initial
	begin
		handle_x = $fopen("data_x.txt");
		handle_y = $fopen("data_y.txt");
		handle_z = $fopen("data_z.txt");
		handle_w = $fopen("data_w.txt");
	end
	initial clk = 1'b1;      
	always #(`clock_period/2) clk = ~clk;
	initial
		begin
				rst_n = 1'b0;
			#(`clock_period*20+1)
				rst_n = 1'b1;
	
				for(i=0;i<400000;i=i+1'b1)
					begin
                                
						#(`clock_period)
						begin
						$fdisplay(handle_x,"%d",x);//将数据保存为txt文件，十六进制格式的数据，可以通过matlab观察相位图
						$fdisplay(handle_y,"%d",y);
						$fdisplay(handle_z,"%d",z);
					 	$fdisplay(handle_w,"%d",w);
						end
                     
					end
			$stop;
		end

endmodule

用python产生超混沌电路效果代码

import numpy as np  
import matplotlib.pyplot as plt  
def lorenz(x):
    """
    计算Lorenz系统的导数。
    
    Lorenz系统是一组用于描述大气对流简化模型的非线性微分方程。这个系统展示了混沌理论的典型行为，即对初始条件的极端敏感性。
    
    参数:
    x (list): 系统变量 [x, y, z, w] 的列表。
    a (float): 系统参数，代表系统的不同特性。
    b (float): 系统参数，通常与系统的耗散率相关。
    c (float): 系统参数，影响系统的非线性行为。
    r (float): 系统参数，代表外部强迫或驱动项。
    
    返回:
    numpy.array: 导数 [dx, dy, dz, dw] 的数组，描述系统如何随时间变化。
    """
    # 计算x方向的导数
    dx = a * (x[1] - x[0]) + x[3]  
    # 计算y方向的导数
    dy = x[0] * (c - x[2]) - x[1]  
    # 计算z方向的导数
    dz = x[0] * x[1] - b * x[2]  
    # 计算w方向的导数，这是原Lorenz系统中没有的，为了扩展或特定应用而添加
    dw = -x[1] * x[2] + r * x[3]  
    # 返回所有方向导数的数组
    return np.array([dx, dy, dz, dw])  
 
print(np.array([0.1, 0, 0, 0]))
import numpy as np

def runge_kutta_2(f, dt, n):
    """
    使用改良欧拉法（二阶龙格-库塔方法）求解微分方程的数值解。
    
    参数:
    f: 微分方程的函数，接受一个数组作为输入，并返回微分方程的导数数组。
    dt: 时间步长，用于控制每一步的大小。
    n: 迭代步数，决定解决问题的总时间步数。
    
    返回:
    数组，包含每一步的结果，每一行代表一个时间步的结果。
    """
    results = [np.array([0.1, 0, 0, 0])]  # 使用列表来存储每一步的结果
    y = np.array([0.1, 0, 0, 0])  # 初始条件
    for i in range(n):
        k1 = dt * f(y)  # 显式欧拉法的一步预测
        y_pred = y + k1  # 预测值
        k2 = dt * f(y_pred)  # 使用预测值计算导数
        y = y + 0.5 * (k1 + k2)  # 更新y值
        results.append(y)  # 存储结果
    return np.array(results)  # 将结果转换为numpy数组并返回

 
# 参数设置  产生随机密码子
a = np.random.uniform(2, 16)  
b = np.random.uniform(2, 4)  
c = np.random.uniform(25, 30)  
r = -0.5 #np.random.uniform(-1, 0)  
print(a, b, c, r)  
  
# 定义时间步长
dt = 0.01  
 
# 定义模拟或计算的步数
num_steps = 10000  
  
 
# 使用改良欧拉法离散化 Lorenz 系统  
x = runge_kutta_2(lambda x: lorenz(x), dt, num_steps)  
  
# 绘制轨迹  
fig = plt.figure()  
ax = fig.add_subplot(121, projection='3d')  
ax.set_title('Lorenz System - x, y, z')  
ax.plot(x[:, 0], x[:, 1], x[:, 2], lw=1)  
ax.set_xlabel('x')  
ax.set_ylabel('y')  
ax.set_zlabel('z')  
  
ax = fig.add_subplot(122, projection='3d')  
ax.set_title('Lorenz System - x, y, w')  
ax.plot(x[:, 0], x[:, 1], x[:, 3], lw=1)  
ax.set_xlabel('x')  
ax.set_ylabel('y')  
ax.set_zlabel('w')  
  
plt.show()

用python检验仿真的数据结果

import numpy as np  
import matplotlib.pyplot as plt  
  
# 加载数据  
def load_data(filename):  
    with open(filename, 'r') as file:  
        data = np.fromfile(file, sep='\n', dtype=int)  
    return data  
  
# 加载所有数据文件  
data_x = load_data('data_x.txt')  
data_y = load_data('data_y.txt')  
data_z = load_data('data_z.txt')  
data_w = load_data('data_w.txt')   
n = len(data_x)  
  
# 创建一个空的numpy数组来存储组合后的数据  
combined_data = np.zeros((n, 4), dtype=int)  
combined_data[:, 0] = data_x  
combined_data[:, 1] = data_y  
combined_data[:, 2] = data_z  
combined_data[:, 3] = data_w  

# 使用 combined_data 来绘图  
fig = plt.figure()  
  
# 第一个3D图：x, y, z  
ax = fig.add_subplot(121, projection='3d')  
ax.set_title('3D Trajectory - x, y, z')  
ax.plot(combined_data[:, 0], combined_data[:, 1], combined_data[:, 2], lw=1)  
ax.set_xlabel('x')  
ax.set_ylabel('y')  
ax.set_zlabel('z')  

# 第二个3D图：x, y, w  
ax = fig.add_subplot(122, projection='3d')  
ax.set_title('3D Trajectory - x, y, w')  
ax.plot(combined_data[:, 0], combined_data[:, 1], combined_data[:, 3], lw=1)  
ax.set_xlabel('x')  
ax.set_ylabel('y')  
ax.set_zlabel('w')  
  
plt.tight_layout()  # 调整子图布局以避免重叠  
plt.show()