基于混沌映射的文本加密算法研究（四）

• 前言
• 一、传统DES密码算法
• 二、典型的文本混沌加密算法
• • 1、Logistic映射
  • 2、Henon映射
  • 3、Lorenz映射
  • 4、Hodgkin-Huxley模型
• 代码

前言

前面的文章介绍了混沌的基础知识，包括三个经典的混沌映射以及Hodgkin-Huxley模型。本文将介绍基于混沌映射的文本加密。

本章先介绍基于传统DES算法的密码系统以及典型的基于混沌映射的文本密码系统，使用前文提出的经典混沌映射和Hodgkin-Huxley模型对文本进行混沌加密，并基于第二章中提出的密码系统性能分析指标对不同算法的安全性进行比较分析。针对典型的文本混沌加密算法不具有密文敏感性的缺陷，设计了一个改进的文本混沌加密算法，并以Hodgkin-Huxley模型为例，对典型的文本混沌加密算法以及改进的文本混沌加密算法进行了比较分析。

一、传统DES密码算法

DES算法是对称分组密码中最为经典的代表，也是NIST建议的文本加密标准算法。DES算法的加解密流程如图5-1所示。

网络上有很多关于DES算法的文章，此处不再介绍。

本文对Txt文本文件进行加密解密，使用的明文文本为马丁\bullet路德\bullet金于1963年发表的著名演讲文本《I Have a Dream》（《我有一个梦想》），文本如图5-3所示。按从左到右、从上到下的顺序将明文的ASCII码展开为一个大小为8866*1的一维列向量，其ASCII码的直方图如图所示：


由图5-4可知，明文文本的ASCII码分布主要集中在32（空格）以及97_{122（小写字母a}z）之间，具有较为明显的非均匀分布特征。
设DES算法的初始密钥为K=133457799BBCDFF1（十六进制），采用图5-1所示的DES算法，在Matlab中进行实验，实验结果如图5-5和图5-6所示：

由图5-5和图5-6可知，传统的DES算法能够达到通过加密隐藏原文信息、通过解密正确恢复原文信息的目的。

二、典型的文本混沌加密算法

基于混沌序列密码的典型混沌文本密码系统如图5-7所示：

在图5-7中，置乱算法采用一维无重复置乱算法，扩散算法采用基于有限域GF\left(257\right)的乘法运算的扩散算法。
下面对基于不同混沌映射的文本混沌加密进行分析。首先使用第三章中介绍的经典混沌映射对文本进行混沌加密，然后使用第四章中介绍的Hodgkin-Huxley模型对文本进行混沌加密。

1、Logistic映射

2、Henon映射

3、Lorenz映射

4、Hodgkin-Huxley模型

代码

由于篇幅原因，仅以HH模型为例子给出代码。

% FIPS 140-1随机特性测试
clear
clc

sigma = 10; alpha = 28;  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);
s = Y1(end-20000+1:end,1);
s = mod(floor((s+1000)*pow2(16)),2);

test = zeros(1,4);
% 单比特测试
k1 = sum(s); k0 = length(s) - k1;
if (9725 <= k1) && (k1 <= 10725) && (9725 <= k0) && (k0 <= 10725)
    test(1) = 1;
end

% 扑克测试
f0 = zeros(1,16);
for i = 1:5000
    v = s(4*i-3)*8 + s(4*i-2)*4 + s(4*i-1)*2 + s(4*i);
    f0(v+1) = f0(v+1) + 1;
end
p = sum(f0.^2)*16/5000 - 5000;
if (p >= 2.16) && (p <= 46.17)
    test(2) = 1;
end

% 游程测试
run0 = zeros(1,100);
run1 = zeros(1,100);
t0 = 0; t1 = 0;
for i = 1:20000
    if i == 1
        if s(i) == 1
            t1 = 1;
        else
            t0 = 1;
        end
    elseif i == 20000
        if s(i) == 1 && t1 > 0
            t1 = t1 + 1;
            run1(t1) = run1(t1) + 1;
        end
        if s(i) == 0 && t0 > 0
            t0 = t0 + 1;
            run0(t0) = run0(t0) + 1;
        end
        if s(i) == 1 && t1 == 0
            run1(1) = run1(1) + 1;
            run0(t0) = run0(t0) + 1;
        end
        if s(i) == 0 && t0 == 0
            run0(1) = run0(1) + 1;
            run1(t1) = run1(t1) + 1;
        end
    else
        if s(i) == 1 && t1 > 0
            t1 = t1 + 1;
        end
        if s(i) == 0 && t0 > 0
            t0 = t0 + 1;
        end
        if s(i) == 1 && t1 == 0
            t1 = 1;
            run0(t0) = run0(t0) + 1;
            t0 = 0;
        end
        if s(i) == 0 && t0 == 0
            t0 = 1;
            run1(t1) = run1(t1) + 1;
            t1 = 0;
        end
    end
end
run0_1 = run0(1); run0_2 = run0(2); run0_3 = run0(3); run0_4 = run0(4);
run0_5 = run0(5); run0_6 = sum(run0(6:100)); run0_26 = sum(run0(26:100));
run1_1 = run1(1); run1_2 = run1(2); run1_3 = run1(3); run1_4 = run1(4);
run1_5 = run1(5); run1_6 = sum(run1(6:100)); run1_26 = sum(run1(26:100));
if (run0_1 >= 2315) && (run0_1 <= 2685) && (run1_1 >= 2315) && (run1_1 <= 2685)
    if (run0_2 >= 1114) && (run0_2 <= 1386) && (run1_2 >= 1114) && (run1_2 <= 1386)
        if (run0_3 >= 527) && (run0_3 <= 723) && (run1_3 >= 527) && (run1_3 <= 723)
            if (run0_4 >= 240) && (run0_4 <= 384) && (run1_4 >= 240) && (run1_4 <= 384)
                if (run0_5 >= 103) && (run0_5 <= 209) && (run1_5 >= 103) && (run1_5 <= 209)
                    if (run0_6 >= 103) && (run0_6 <= 207) && (run1_6 >= 103) && (run1_6 <= 207)
                        test(3) = 1;
                    end
                end
            end
        end
    end
end

if (run0_26 == 0) && (run1_26 == 0)
    test(4) = 1;
end

clearvars -except k1 p run0_* run1_* test

%%
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 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

%% HH_Lorenz混沌加密系统 性能分析
clear
clc
% 读取明文
plaintext = importdata('D:\学习资料\本科毕业论文\代码\加密\plain.txt');
[m,~] = size(plaintext);
ascsize = zeros(m,1);
plaintext_asc = [];
for i = 1:m
    text = char(plaintext(i));
    asc = abs(text);
    ascsize(i) = length(asc);
    if i > 1
        ascsize(i) = ascsize(i) + ascsize(i-1);
    end
    plaintext_asc = [plaintext_asc asc];
end
plaintext_asc = plaintext_asc';

clearvars -except plaintext* ascsize m

%%
% 密钥K=[alpha,V0,m0,n0,h0]
K = [28, 0.00002, 0.05293, 0.31768, 0.59612];
% 加密
tic
ciphetext_asc = Encrypt(plaintext_asc,K);
toc
ciphetext = cell(m,1);
for i = 1:m
    if i == 1
        asc = ciphetext_asc(1:ascsize(1))';
    else
        asc = ciphetext_asc(ascsize(i-1)+1:ascsize(i))';
    end
    text = char(asc);
    ciphetext(i) = cellstr(text);
end
ciphetext = char(ciphetext);

% 解密
tic
plaintext2_asc = Decrypt(ciphetext_asc,K);
toc
plaintext2 = cell(m,1);
for i = 1:m
    if i == 1
        asc = plaintext2_asc(1:ascsize(1))';
    else
        asc = plaintext2_asc(ascsize(i-1)+1:ascsize(i))';
    end
    text = char(asc);
    plaintext2(i) = cellstr(text);
end
plaintext2 = char(plaintext2);

clearvars -except K plaintext* ciphetext* ascsize m

%% 加密/解密速度
K = [28, 0.00002, 0.05293, 0.31768, 0.59612];
N = 100;
time1 = zeros(N,1); time2 = zeros(N,1);
for j = 1:N
    time1(j) = cputime;
    ciphe = Encrypt(plaintext_asc,K);
    time1(j) = cputime - time1(j);
    time2(j) = cputime;
    Decrypt(ciphe,K);
    time2(j) = cputime - time2(j);
end
En_time = sum(time1)/N;
De_time = sum(time2)/N;

clearvars j N time*

%% 直方图卡方检验
figure
hist(plaintext_asc,256);
figure
hist(ciphetext_asc,256);
figure
hist(plaintext2_asc,256);

fp = hist(plaintext_asc,256);
fc = hist(ciphetext_asc,256);
g = ascsize(end)/256;
% 卡方统计量的值
chai1 = sum((fp-g).^2)/g;
chai2 = sum((fc-g).^2)/g;

clearvars M N fp fc g


%% NPCR和UACI值
NU = NPCRUACI(plaintext_asc, ciphetext_asc);

%%
% 给定文本和随机文本的期望值
NUExpect = NPCRUACIExpect(plaintext_asc);
% 两个随机文本的理论值
NUTheory = NPCRUACITheory();

%% 加密系统的密钥敏感性
N = 100;
alpha = roundn(rand(1,N)*24, -2) + 66;
V0 = roundn(rand(1,N)*10, -10);
m0 = roundn(rand(1,N)*0.1, -10);
n0 = roundn(rand(1,N)*0.1, -10) + 0.3;
h0 = roundn(rand(1,N)*0.2, -10) + 0.4;
NU_alpha = zeros(1,2);
NU_V0 = zeros(1,2);
NU_m0 = zeros(1,2);
NU_n0 = zeros(1,2);
NU_h0 = zeros(1,2);
for i = 1:N
    K1 = [alpha(i), V0(i), m0(i), n0(i), h0(i)];
    ciphe_asc1 = Encrypt(plaintext_asc, K1);
    K2 = [alpha(i)+10^(-10), V0(i), m0(i), n0(i), h0(i)];
    if K2(1) >= 90
        K2(1) = alpha(i) - 10^(-10);
    end
    ciphe_asc2 = Encrypt(plaintext_asc, K2);    
    K3 = [alpha(i), V0(i)+10^(-10), m0(i), n0(i), h0(i)];
    if K3(2) >= 10
        K3(2) = V0(i) - 10^(-10);
    end
    ciphe_asc3 = Encrypt(plaintext_asc, K3);    
    K4 = [alpha(i), V0(i), m0(i)+10^(-10), n0(i), h0(i)];
    if K4(3) >= 0.1
        K4(3) = m0(i) - 10^(-10);
    end
    ciphe_asc4 = Encrypt(plaintext_asc, K4);        
    K5 = [alpha(i), V0(i), m0(i), n0(i)+10^(-10), h0(i)];
    if K5(4) >= 0.4
        K5(4) = n0(i) - 10^(-10);
    end
    ciphe_asc5 = Encrypt(plaintext_asc, K5);       
    K6 = [alpha(i), V0(i), m0(i), n0(i), h0(i)+10^(-10)];
    if K6(5) >= 0.6
        K6(5) = h0(i) - 10^(-10);
    end
    ciphe_asc6 = Encrypt(plaintext_asc, K6);       
    
    NU_alpha = NU_alpha + NPCRUACI(ciphe_asc1, ciphe_asc2);
    NU_V0 = NU_V0 + NPCRUACI(ciphe_asc1, ciphe_asc3);
    NU_m0 = NU_m0 + NPCRUACI(ciphe_asc1, ciphe_asc4);
    NU_n0 = NU_n0 + NPCRUACI(ciphe_asc1, ciphe_asc5);
    NU_h0 = NU_h0 + NPCRUACI(ciphe_asc1, ciphe_asc6);
end
NU_alpha = NU_alpha/N;
NU_V0 = NU_V0/N;
NU_m0 = NU_m0/N;
NU_n0 = NU_n0/N;
NU_h0 = NU_h0/N;

%% 明文的敏感性分析
N = 50;
alpha = roundn(rand(1,N)*24, -2) + 66;
V0 = roundn(rand(1,N)*10, -10);
m0 = roundn(rand(1,N)*0.1, -10);
n0 = roundn(rand(1,N)*0.1, -10) + 0.3;
h0 = roundn(rand(1,N)*0.2, -10) + 0.4;
NU = zeros(N,2);
NU_plain = zeros(1,2);
for i = 1:N
    K1 = [alpha(i), V0(i), m0(i), n0(i), h0(i)];
    ciphe_asc1 = Encrypt(plaintext_asc, K1);
    ixy = mod(floor(rand*10^10), ascsize(end))+1;
    plain_asc = plaintext_asc;
    plain_asc(ixy) = mod(plain_asc(ixy)+1,256);
    ciphe_asc2 = Encrypt(plain_asc, K1);
    NU(i,:) = NPCRUACI(ciphe_asc1,ciphe_asc2);
    NU_plain = NU_plain + NPCRUACI(ciphe_asc1, ciphe_asc2);
end
NU_plain = NU_plain/N;

clearvars i* K* ciphe_asc* plain_asc N

%% 密文的敏感性分析
N = 100;
alpha = roundn(rand(1,N)*24, -2) + 66;
V0 = roundn(rand(1,N)*10, -10);
m0 = roundn(rand(1,N)*0.1, -10);
n0 = roundn(rand(1,N)*0.1, -10) + 0.3;
h0 = roundn(rand(1,N)*0.2, -10) + 0.4;
NU = zeros(N,2);
NU_ciphe = zeros(1,2);
for i = 1:N
    K1 = [alpha(i), V0(i), m0(i), n0(i), h0(i)];
    ciphe_asc = Encrypt(plaintext_asc, K1);
    ixy = mod(floor(rand*10^10), ascsize(end))+1;
    ciphe_asc(ixy) = mod(ciphe_asc(ixy)+1,256);
    plain_asc = Decrypt(ciphe_asc, K1);
    NU(i,:) = NPCRUACI(plaintext_asc,plain_asc);
    NU_ciphe = NU_ciphe + NPCRUACI(plaintext_asc, plain_asc);
end
NU_ciphe = NU_ciphe/N;

clearvars i* ciphe_asc plain_asc K1 N

%% 信息熵
alpha = roundn(rand*24, -2) + 66;
V0 = roundn(rand*10, -10);
m0 = roundn(rand*0.1, -10);
n0 = roundn(rand*0.1, -10) + 0.3;
h0 = roundn(rand*0.2, -10) + 0.4;
ciphe_asc = Encrypt(plaintext_asc, [alpha, V0, m0, n0, h0]);
Hplain = Entropy(plaintext_asc);
Hciphe = Entropy(ciphe_asc);

clearvars ciphe_asc

% HH混沌加密系统加密函数Encrypt
function ciphetext_asc = Encrypt(P,K)
% P为明文，K为密钥
num1 = length(P);
num2 = 3*num1;
s = zeros(1,num2);

% 生成混沌序列
sigma = 10; alpha = K(1); beta= 8/3;
K0 = [K(2), K(3), K(4), K(5), 1, 1, 1];
[T,Y] = ode45(@(t,x)HHfun(t,x,sigma,alpha,beta),[0,5000],K0);
s = Y(end-num2+1:end,1);

% 置乱算法密码流
X = mod(floor(s(1:num1)*10^10),num1)+1;
[~,idx] = unique(X);
X1 = zeros(1,num1);
X1(1:length(idx)) = X(sort(idx));
X1(length(idx)+1:num1) = setdiff(1:num1,X1);
X = X1;

% 扩散算法密码流
S = mod(floor(s(num1+1:3*num1)*pow2(16)),256);
S1 = S(1:num1); S2 = S(num1+1:2*num1);

% 生成乘法运算表
TBL = GF257Table();

% 置乱算法
A = P;
for i = 1:floor(num1/2)
    t = A(X(i));
    A(X(i)) = A(X(num1-i+1));
    A(X(num1-i+1)) = t;
end

% 扩散算法
B = zeros(num1,1); C = zeros(num1,1);
B0 = 0;
B(1) = LookUpGF257(B0,S1(1),A(1),TBL);
for i = 2:num1
    B(i) = LookUpGF257(B(i-1),S1(i),A(i),TBL);
end
C0 = 0;
C(num1) = LookUpGF257(C0,S2(num1),B(num1),TBL);
for i = num1-1:-1:1
    C(i) = LookUpGF257(C(i+1),S2(i),B(i),TBL);
end
ciphetext_asc = C;
end

% Logistic混沌加密系统解密函数Decrypt
function plaintext2_asc = Decrypt(C,K)
% C为密文，K为密钥

num1 = length(C);
num2 = 3*num1;
s = zeros(1,num2);

% 生成混沌序列
sigma = 10; alpha = K(1); beta= 8/3;
K0 = [K(2), K(3), K(4), K(5), 1, 1, 1];
[T,Y] = ode45(@(t,x)HHfun(t,x,sigma,alpha,beta),[0,5000],K0);
s = Y(end-num2+1:end,1);

% 置乱算法密码流
X = mod(floor(s(1:num1)*10^10),num1)+1;
[~,idx] = unique(X);
X1 = zeros(1,num1);
X1(1:length(idx)) = X(sort(idx));
X1(length(idx)+1:num1) = setdiff(1:num1,X1);
X = X1;

% 扩散算法密码流
S = mod(floor(s(num1+1:3*num1)*pow2(16)),256);
S1 = S(1:num1); S2 = S(num1+1:2*num1);

% 生成除法运算表
[TBL1,TBL2] = GF257TableEx();

% 扩散算法的逆算法
A = C;
D = zeros(num1,1); E = zeros(num1,1);
A0 = 0;
D(num1) = LookUpGF257Ex2(A(num1),A0,S2(num1),TBL1,TBL2);
for i = num1-1:-1:1
    D(i) = LookUpGF257Ex2(A(i),A(i+1),S2(i),TBL1,TBL2);
end
E0 = 0; E(1) = LookUpGF257Ex2(D(1),E0,S1(1),TBL1,TBL2);
for i = 2:num1
    E(i) = LookUpGF257Ex2(D(i),D(i-1),S1(i),TBL1,TBL2);
end

% 置乱算法的逆算法
for i = 1:floor(num1/2)
    t = E(X(i));
    E(X(i)) = E(X(num1-i+1));
    E(X(num1-i+1)) = t;
end
plaintext2_asc = E;
end

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

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

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

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

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

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

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

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

function gvalue = g(x1,x2,x3,alpha)
gvalue = alpha*x1 - x1*x3 - x2;
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

% 基于GF(257)域的乘法运算表
function T = GF257Table()
T = mod(transpose(0:256).*(0:256),257);
end

% 基于GF(257)域的乘法和除法运算表
function [T1,T2] = GF257TableEx()
T1 = mod(transpose(0:256).*(0:256),257);
T2 = zeros(257,257);
T2(:,1) = transpose(0:256);
for j = 1:256
    for k = 1:256
        T2(mod(j*k,257)+1,j+1) = k;
    end
end
end

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 lvalue = l(x1,x2,x3,beta)
lvalue = x1*x2 - beta*x3;
end

% 基于GF(257)域的乘法运算函数
function y = LookUpGF257(x0,x1,x2,TBL)
y = TBL(TBL(x0+2,x1+2)+1,x2+2)-1;
end

% 基于GF(257)域的除法运算函数
function y = LookUpGF257Ex2(x0,x1,x2,TBL1,TBL2)
t = TBL1(x1+2,x2+2)+1;
y = TBL2(x0+2,t)-1;
end

% 计算NPCR、UACI和BACI值
function u = NPCRUACI(A,B)
m = length(A);
u = zeros(1,2);

u(1) = sum(sum(A ~= B))/m*100;
u(2) = sum(sum(abs(A - B)))/(255*m)*100;
end

% 计算给定文本你与随机文本的NPCR期望值和UACI期望值
function u = NPCRUACIExpect(A)
m = length(A);
u = zeros(1,2);

u(1) = 255/256*100;

tot_n = zeros(1,256);
tot_s = 0:255;
for i = 1:m
    for k = 0:255
        if k <= A(i)
            tot_n(k+1) = tot_n(k+1)+1;
        end
        if k <= 255 - A(i)
            tot_n(k+1) = tot_n(k+1)+1;
        end
    end
end
tot_n(1) = tot_n(1)/2;
u(2) = sum(tot_s.*tot_n)/sum(tot_n)/255*100;
end

% 两个随机文本的NPCR和UACI理论值
function u = NPCRUACITheory()
u = zeros(1,2);
u(1) = 255/256*100;
a = 255:-1:1; b = 1:255;
u(2) = sum(a.*b)*2/(256^2*255)*100;
end