基于混沌系统的文本加密算法研究系列

最新推荐文章于 2023-04-11 21:19:01 发布

单身苟阿富汗

最新推荐文章于 2023-04-11 21:19:01 发布

阅读量1k

点赞数 2

文章标签：算法 matlab 加密解密

本文链接：https://blog.csdn.net/qq_42665419/article/details/120402498

版权

基于混沌系统的文本加密算法研究（三）

前言
一、 Hodgkin-Huxley模型的数学形式
- 2、Hodgkin-Huxley模型的混沌分析
- - （1）外部电流为混沌序列
  - （2）外部电流为周期性刺激电流
总结
代码
- 1、外部电流为混沌序列
- 2、外部电流为周期性刺激电流

前言

前面的文章介绍了混沌的基础知识以及三个经典的混沌映射，本文将介绍著名的Hodgkin-Huxley模型。
（传送门：基于混沌系统的文本加密算法研究（一）——混沌及混沌加密的基础知识）
（传送门：基于混沌系统的文本加密算法研究（二）——经典混沌映射）

1952年英国生理学家Alan Hodgkin和Andrew Huxley在对乌贼神经轴突电生理特性研究的基础上，提出了Hodgkin-Huxley（HH）模型。本章将介绍经典的HH模型以，并对HH模型的混沌性质进行研究。

一、 Hodgkin-Huxley模型的数学形式

20世纪50年代，Hodgkin和Huxley利用电压钳位技术对乌贼巨轴突的生理特性进行了深入的研究，在此基础上于1952年总结提出了经典的HH模型。该模型随后成为了很多不同生理结构中神经细胞的模型的基础。经典HH模型由四个微分方程组成，其表达式如式(4-1)所示：
$\left\{ \begin{array}{l}\frac{dV}{dt}=\frac{1}{C_M}[I-\bar{g}_{Na}m^3h(V-V_{Na})-\bar{g}_Kn^4(V-V_K)-\bar{g}_l(V-V_l)]\\ \frac{dm}{dt}=\alpha_m(1-m)-\beta_mm\\ \frac{dh}{dt}=\alpha_h(1-h)-\beta_hh\\ \frac{dn}{dt}=\alpha_n(1-n)-\beta_nn\\ \end{array} \right. \tag{1}$
其中，V为跨膜电压；m、h、n分别表示钠离子激活变量、钠离子抑制变量和钾离子激活变量；I为外刺激电流； $V_{Na}、V_{K}、V_{l}$ 分别表示钠离子平衡电位、钾离子平衡电位和漏离子平衡电位； $\bar{g}_{Na}、\bar{g}_{K}、\bar{g}_{l}$ 分别表示钠离子通道最大电导、钾离子通道最大电导和漏离子通道最大电导； $\alpha_m、\beta_m、\alpha_h、\beta_h、\alpha_n和\beta_n$ 的表达式由实验数据推导得出，如式（2）所示。
$\alpha_m(V)=\frac{\bar{\alpha}_m(V-\bar{V}_m)}{1-exp[-(V-\bar{V}_m)/K_{\alpha_m}]}$ $\beta_m(V)=\bar{\beta}_mexp(-V/K_{\beta_m})$ $\alpha_h(V)=\bar{\alpha}_hexp(-V/K_{\alpha_h}) \tag{2}$ $\beta_h(V)=\frac{\bar{\beta}_h}{1+exp[-(V-\bar{V}_h)/K_{\beta_h}]}$ $\alpha_n(V)=\frac{\bar{\alpha}_n(V-\bar{V}_n)}{1-exp[-(V-\bar{V}_n)/K_{\alpha_n}]}$ $\beta_n(V)=\bar{\beta}_nexp(-V/K_{\beta_n})$ 公式（1）和（2）是温度为 $6.3^{\circ}$ C时的经典HH模型的表达式，其中各参数的取值如下表所示。
在这里插入图片描述

2、Hodgkin-Huxley模型的混沌分析

对HH模型的混沌分析可以从两个方面进行研究。一方面，可以研究外部混沌信号对HH模型的影响，即将混沌系统生成的混沌序列作为电流输入HH模型中；另一方面，可以通过引入周期性电流刺激，对HH模型的混沌性质进行研究。

（1）外部电流为混沌序列

首先对外部电流为混沌序列的情形进行研究。
在这里插入图片描述

（2）外部电流为周期性刺激电流

下面对引入周期性电流刺激的情形进行研究。

在这里插入图片描述

总结

本文介绍了Hodgkin-Huxley模型，对HH模型的的混沌性质进行了研究。通过将由Lorenz系统生成的混沌信号作为外部刺激电流，对Lorenz系统控制参数与神经纤维的放电模式的关系进行分析。通过实验仿真可知，随着控制参数\alpha的增加，神经纤维的放电模式经历了从静息状态到混沌放电的转化过程；且当Lorenz系统处于混沌状态时，以Lorenz系统生成的混沌信号作为刺激电流，HH模型也将进入混沌状态。

代码

1、外部电流为混沌序列

%% ode45
clear
clc
sigma = 10; alpha = 92;  beta= 8/3;
K0 = [0.00002, 0.05293, 0.31768, 0.59612, 1, 1, 1];
[T1,Y1] = ode45(@(t,x)HHfun(t,x,sigma,alpha,beta),[0,5000],K0);
plot(T1,Y1(:,1),'k')
xlim([4900,5000])
xlabel('time/ms')
ylabel('V')
figure
plot(T1,Y1(:,5),'k')
xlim([4900,5000])
xlabel('time/ms')
ylabel('x')

%% 极值点图(ode45)
clear
sigma = 10; alpha = 0:0.1:100;  beta= 8/3;
K0 = [0.00002, 0.05293, 0.31768, 0.59612, 1, 1, 1];
point = [];
for j = 1:length(alpha)
    [T,Y1] = ode45(@(t,x)HHfun(t,x,sigma,alpha(j),beta),[0,2000],K0);
    V = Y1(T>1900,1);
    maxmin1 = findpeaks(V);
    maxpoint = max(V);
    maxmin2 = findpeaks(maxpoint - V);
    for k = 1:length(maxmin1)
        point = [point maxmin1(k)*i+alpha(j)];
    end
    for k = 1:length(maxmin2)
        point = [point (maxpoint - maxmin2(k))*i+alpha(j)];
    end
end
plot(point,'k.')
xlabel('\alpha')
ylabel('极值点/mV')

%% 最大Lyapunov指数(ode45)
clear
clc
sigma = 10; alpha = 20:0.01:100;  beta= 8/3;
d0 = 10^(-10);
N0 = 2000;
N = 1000;
Z = [];
for k = 1:length(alpha)
    Lya = 0;
    K0 = [0.00002, 0.05293, 0.31768, 0.59612, 1, 1, 1];
    K1 = [0.00002, 0.05293+10^(-10), 0.31768, 0.59612, 1 , 1, 1];
    for j = 1:N0+N
        [T1,Y1] = ode45(@(t,x)HHfun(t,x,sigma,alpha(k),beta),[0,1],K0);
        
        [T2,Y2] = ode45(@(t,x)HHfun(t,x,sigma,alpha(k),beta),[0,1],K1);
        d1 = sqrt((Y1(end,1) - Y2(end,1))^2 + (Y1(end,2) - Y2(end,2))^2 +(Y1(end,3) - Y2(end,3))^2 +(Y1(end,4) - Y2(end,4))^2 +...
            (Y1(end,5) - Y2(end,5))^2 +(Y1(end,6) - Y2(end,6))^2 + (Y1(end,7) - Y2(end,7))^2);
        K0 = Y1(end,:);
        for i = 1:7
            K1(i) = Y1(end,i) + (d0/d1)*(Y2(end,i) - Y1(end,i));
        end
        if j > N0
            Lya = Lya + log(d1/d0);
        end
    end
    Z = [Z Lya/N];
end
plot(alpha,Z,'k')
hold on
xlabel('\alpha')
ylabel('最大Lyapunov指数')
plot(alpha,zeros(1,length(alpha)),'r-.')

%%
function df = HHfun(t,x,sigma,alpha,beta)
df = zeros(7,1);
df(1) = x(5) - (120*x(2)^3*x(4)*(x(1) - 115) + 36*x(3)^4*(x(1) + 12) + 0.3*(x(1) - 10.599));
df(2) = alpham(x(1))*(1-x(2)) - betam(x(1))*x(2);
df(3) = alphan(x(1))*(1-x(3)) - betan(x(1))*x(3);
df(4) = alphah(x(1))*(1-x(4)) - betah(x(1))*x(4);
df(5) = sigma*(x(6) - x(5));
df(6) = alpha*x(5) - x(5)*x(7) - x(6);
df(7) = x(5)*x(6) - beta*x(7);
end


function fvalue = f(x1,x2,sigma)
fvalue = sigma*(x2 - x1);
end

function gvalue = g(x1,x2,x3,alpha)
gvalue = alpha*x1 - x1*x3 - x2;
end

function lvalue = l(x1,x2,x3,beta)
lvalue = x1*x2 - beta*x3;
end

function fvalue = f1(m,n,h,V,x1)
fvalue = x1 - (120*m^3*h*(V-115) + 36*n^4*(V+12) + 0.3*(V-10.599));
end

function gvalue = g1(m,V)
gvalue = alpham(V)*(1-m) - betam(V)*m;
end

function gvalue = g2(n,V)
gvalue = alphan(V)*(1-n) - betan(V)*n;
end

function gvalue = g3(h,V)
gvalue = alphah(V)*(1-h) - betah(V)*h;
end

function am = alpham(V)
am = 0.1*(V - 25)/(1 - exp((25 - V)/10));
end

function bm = betam(V)
bm = 4*exp(-V/18);
end

function an = alphan(V)
an = 0.01*(V - 10)/(1 - exp((10 - V)/10));
end

function bn = betan(V)
bn = 0.125*exp(-V/80);
end

function ah = alphah(V)
ah = 0.07*exp(-V/20);
end

function bh = betah(V)
bh = 1/(1 + exp((30 - V)/10));
end

2、外部电流为周期性刺激电流

%% HH模型
clear
clc
A = roundn(rand*2, -2) + 4;
f = roundn(rand*20 , -1) + 120;
V0 = roundn(rand*110, -10) - 10;
m0 = roundn(rand*1, -10);
n0 = roundn(rand*0.45, -10) + 0.3;
h0 = roundn(rand*0.5, -10) + 0.1;
x0 = [V0, m0, n0, h0, 0];
[T,Y] = ode45(@(t,x)HHfun(t,x,A,f),[0,500],x0);
V = Y(:,1);
m = Y(:,2);
n = Y(:,3);
h = Y(:,4);
plot(T,V)

%% HH模型（f）
clear
clc
x0 = [0.00002, 0.05293, 0.31768, 0.59612, 0];
A = 5; f = 130;
[T,Y] = ode45(@(t,x)HHfun(t,x,A,f),[0,500],x0);
V = Y(:,1);
m = Y(:,2);
n = Y(:,3);
h = Y(:,4);
% 相空间图
figure 
plot(V,m) xlabel('V') ylabel('m')
figure 
plot(V,n) xlabel('V') ylabel('n')
figure 
plot(V,h) xlabel('V') ylabel('h')
% plot(T,V)
% I = A*(cos(2*pi*f*T/1000) + 1)/2 - 30;
% plot(T,V,'k')
% hold on
% plot(T,I,'k')
% set(gca,'ytick',[-30,0,50,100])
% set(gca,'yticklabel',[0,0,50,100])
% xlabel('time/ms')
% ylabel('V/mV')
% grid on
% xlim([-10,20])


%% 极值点分布(f)
clear
clc
A = 5; f = 0:0.1:150;
x0 = [0.00002, 0.05293, 0.31768, 0.59612, 0];
point = [];
for j = 1:length(f)
    [T,Y] = ode45(@(t,x)HHfun(t,x,A,f(j)),[0,500],x0);
    V = Y(:,1);
    maxmin1 = findpeaks(V);
    maxpoint = max(V);
    maxmin2 = findpeaks(maxpoint - V);
    clearvars i
    for k = 1:length(maxmin1)
        point = [point maxmin1(k)*i+f(j)];
    end
    for k = 1:length(maxmin2)
        point = [point (maxpoint - maxmin2(k))*i+f(j)];
    end
end
plot(point,'k.')
set(gca,'xtick',0:10:150)
xlabel('f/Hz')
ylabel('V/mV')

%% 最大Lyapunov指数图(f)
clear
clc
A = 5; f = 120:0.01:145;
t = 500;
N = 2000;
d0 = 1e-10;
Z = [];
for j = 1:length(f)
    Lya = 0;
    x0 = [0.00002, 0.05293, 0.31768, 0.59612, 0];
    x1 = x0;
    x1(1) = x1(1) + d0;
    for k = 1:N
        [T1,Y1] = ode45(@(t,x)HHfun(t,x,A,f(j)),[0,1],x0);
        [T2,Y2] = ode45(@(t,x)HHfun(t,x,A,f(j)),[0,1],x1);
        x0 = Y1(end,:); x1 = Y2(end,:);
        d1 = sqrt((x0(1) - x1(1))^2 + (x0(2) - x1(2))^2 + (x0(3) - x1(3))^2 + (x0(4) - x1(4))^2 + (x0(5) - x1(5))^2);
        x1 = x0 + (d0/d1)*(x1 - x0);
        if k > t
            Lya = Lya + log(d1/d0);
        end
    end
    Z = [Z Lya/(k-t)];
end
plot(f,Z,'k')
hold on
plot(f,zeros(1,length(f)),'r-')
xlabel('f/Hz')
ylabel('Lyapunov')
% xlim([122,146])


%% HH模型（A）
clear
clc
x0 = [0.00002, 0.05293, 0.31768, 0.59612, 0];
A = 5; f = 130;
[T,Y] = ode45(@(t,x)HHfun(t,x,A,f),[0,1000],x0);
V = Y(:,1);
V = Y(:,1);
m = Y(:,2);
n = Y(:,3);
h = Y(:,4);
% 相空间图
figure 
plot(V,m)
xlabel('V')
ylabel('m')
figure 
plot(V,n)
xlabel('V')
ylabel('n')
figure 
plot(V,h)
xlabel('V')
ylabel('h')
% 
% plot(T,V)
% I = A*(cos(2*pi*f*T/1000) + 1)/2 - 30;
% plot(T,V,'k')
% hold on
% plot(T,I,'k')
% % set(gca,'ytick',[-6,-2,-1,0,1,2,3,4])
% % set(gca,'yticklabel',[0,-2,-1,0,1,2,3,4])
% set(gca,'ytick',[-30,0,50,100])
% set(gca,'yticklabel',[0,0,50,100])
% xlabel('time/ms')
% ylabel('V/mV')
% grid on

%% 极值点分布（A）
clear
clc
A = 3:0.01:7; f = 130;
x0 = [0.00002, 0.05293, 0.31768, 0.59612, 0];
point = [];
for j = 1:length(A)
    [T,Y] = ode45(@(t,x)HHfun(t,x,A(j),f),[0,1000],x0);
    V = Y(:,1);
    maxmin1 = findpeaks(V);
    maxpoint = max(V);
    maxmin2 = findpeaks(maxpoint - V);
    clearvars i
    for k = 1:length(maxmin1)
        point = [point maxmin1(k)*i+A(j)];
    end
    for k = 1:length(maxmin2)
        point = [point (maxpoint - maxmin2(k))*i+A(j)];
    end
end
plot(point,'k.')
xlabel('A/\muA/cm^{2}')
ylabel('V/mV')

%% 最大Lyapunov指数图（A）
clear
clc
A = 4:0.001:6; f = 130;
t = 500;
N = 1000;
d0 = 1e-10;
Z = [];
for j = 1:length(A)
    Lya = 0;
    x0 = [0.00002, 0.05293, 0.31768, 0.59612, 0];
    x1 = x0;
    x1(1) = x1(1) + d0;
    for k = 1:N
        [T1,Y1] = ode45(@(t,x)HHfun(t,x,A(j),f),[0,1],x0);
        [T2,Y2] = ode45(@(t,x)HHfun(t,x,A(j),f),[0,1],x1);
        x0 = Y1(end,:); x1 = Y2(end,:);
        d1 = sqrt((x0(1) - x1(1))^2 + (x0(2) - x1(2))^2 + (x0(3) - x1(3))^2 + (x0(4) - x1(4))^2 + (x0(5) - x1(5))^2);
        x1 = x0 + (d0/d1)*(x1 - x0);
        if k > t
            Lya = Lya + log(d1/d0);
        end
    end
    Z = [Z Lya/(k-t)];
end
plot(A,Z,'k')
hold on
plot(A,zeros(1,length(A)),'r-')
xlabel('A/\muA/cm^{2}')
ylabel('Lyapunov')
% xlim([122,146])

%% 最大Lyapunov指数图（A、f）
clear
clc
A = 4:0.001:6; f = 130;
t = 500;
N = 1000;
d0 = 1e-10;
Z = [];
for j = 1:length(A)
    Lya = 0;
    x0 = [0.00002, 0.05293, 0.31768, 0.59612, 0];
    x1 = x0;
    x1(1) = x1(1) + d0;
    for k = 1:N
        [T1,Y1] = ode45(@(t,x)HHfun(t,x,A(j),f),[0,1],x0);
        [T2,Y2] = ode45(@(t,x)HHfun(t,x,A(j),f),[0,1],x1);
        x0 = Y1(end,:); x1 = Y2(end,:);
        d1 = sqrt((x0(1) - x1(1))^2 + (x0(2) - x1(2))^2 + (x0(3) - x1(3))^2 + (x0(4) - x1(4))^2 + (x0(5) - x1(5))^2);
        x1 = x0 + (d0/d1)*(x1 - x0);
        if k > t
            Lya = Lya + log(d1/d0);
        end
    end
    Z = [Z Lya/(k-t)];
end
plot(A,Z,'k')
hold on
plot(A,zeros(1,length(A)),'r-')
xlabel('A/\muA/cm^{2}')
ylabel('Lyapunov')
% xlim([122,146])

%%
function df = HHfun(t,x,A,f)
df = zeros(5,1);
df(1) = A*(cos(2*pi*f*x(5)/1000)+1)/2 -(120*x(2)^3*x(4)*(x(1)-115) + 36*x(3)^4*(x(1)+12) + 0.3*(x(1)-10.599));
df(2) = alpham(x(1))*(1-x(2)) - betam(x(1))*x(2);
df(3) = alphan(x(1))*(1-x(3)) - betan(x(1))*x(3);
df(4) = alphah(x(1))*(1-x(4)) - betah(x(1))*x(4);
df(5) = 1;
end

function fvalue = f1(m,n,h,z,V,A,f)
fvalue = A*(cos(2*pi*f*z/1000)+1)/2 -(120*m^3*h*(V-115) + 36*n^4*(V+12) + 0.3*(V-10.599));
end

function gvalue = g1(m,V)
gvalue = alpham(V)*(1-m) - betam(V)*m;
end

function gvalue = g2(n,V)
gvalue = alphan(V)*(1-n) - betan(V)*n;
end

function gvalue = g3(h,V)
gvalue = alphah(V)*(1-h) - betah(V)*h;
end

function am = alpham(V)
am = 0.1*(V - 25)/(1 - exp((25 - V)/10));
end

function bm = betam(V)
bm = 4*exp(-V/18);
end

function an = alphan(V)
an = 0.01*(V - 10)/(1 - exp((10 - V)/10));
end

function bn = betan(V)
bn = 0.125*exp(-V/80);
end

function ah = alphah(V)
ah = 0.07*exp(-V/20);
end

function bh = betah(V)
bh = 1/(1 + exp((30 - V)/10));
end