图像加密这门手艺，做到后期拼的不是花哨的算法数量，而是“你到底用什么手段把像素能量打散”。我最近用Matlab实现了一套基于分数阶傅立叶变换和曲线锯变换的图像加密方案，跑完256×256标准测试图之后，可以说这套组合在统计特性和密钥敏感性上确实比单靠混沌置乱更干净。如果你也在找一种“既能打乱空域位置、又能在变换域做二次保护”的思路，这篇内容应该适合你。
文章里不会只放一个成品代码，我会把为什么选FrFT、曲线锯变换怎么构造、扩散层怎么和复数域配合、解密有哪些参数同步坑，全部拆开讲。适合正在写图像安全相关大作业、课程设计的同学，也适合做论文复现的同行。你不需要事先精通分数阶傅立叶变换，只需要会基本Matlab矩阵操作和一点点信号处理直觉。
1. 为什么我不再用单一的混沌置乱做图像加密
1.1 老方案的问题：置乱能打乱位置，打不散频域能量
先说结论：单纯靠混沌序列打乱像素位置，也就是通常说的“置乱”，确实能让整张图看起来像雪花点，但它在统计特性上有一个隐蔽的软肋——像素值的直方图几乎不变。因为置乱只是重新排列像素坐标，像素数值本身没有改变，所以直方图统计、灰度分布、能量分布这些信息会原封不动地泄露给攻击者。
你可以做一个小实验：把Lena图用Logistic映射生成的索引随机打乱，输出图肉眼完全看不出原图，但一旦对密文做直方图统计，就会看到明显的人脸灰度分布轮廓。这在实际安全性评估里是很掉分的一项，因为直方图泄露意味着明文的统计特征没有被打散。
另一个问题是频域能量。普通置乱在空域里折腾，但图像的相邻像素相关性虽然在空域被破坏了，经过傅立叶变换之后，低频与高频的分布模式仍可能保留原始结构的影子。攻击者用频谱分析或已知明文攻击，很容易反推置乱轨迹的规律。
所以后来做图像加密的人都倾向于“置乱+变换+扩散”三件套。置乱负责打乱位置，变换负责把像素值踢到另一个域去，扩散负责让明文的任何一点微小变化扩散到整张密文。我这次用分数阶傅立叶变换（FrFT）做变换层，用曲线锯变换做置乱层，思路就是这个框架。
1.2 引入FrFT：一个阶数就是一把额外的钥匙
分数阶傅立叶变换在我眼里最迷人的地方是：它把普通傅立叶变换扩展到“分数阶”，你可以理解成在时域和频域之间有一个连续的旋转角度。阶数a=0时信号原样输出，a=1时就是标准傅立叶变换，a=0.5时介于时域和频域之间。对图像做这种变换后，输出是复数矩阵，实部和虚部都携带能量，视觉上完全看不出原始结构。
这个性质对加密来说太有价值了。阶数a不再固定为1，它本身就可以作为一把密钥。解密时必须使用完全相同的阶数做逆变换，哪怕阶数偏差0.0001，解出来的图像就是一团噪声。相比传统固定傅立叶变换，等于白送了一个连续型密钥参数。
另一个好处是FrFT的核函数与阶数绑定，不同的阶数对应不同的能量分布方式。攻击者哪怕拿到了变换后的矩阵，如果不知道阶数和具体离散化方式，也没法把频域数据正确转回空域。这在已知明文攻击模型下比固定变换更抗打。
1.3 曲线锯变换在整套流程里负责什么
标题里的“曲线锯变换”，实现上我采用的是基于空间填充曲线的扫描置乱：先按Zigzag曲线、希尔伯特曲线这类路径，把图像矩阵读成一维序列，再用混沌序列生成一个随机重排索引，把一维序列打乱，最后重排成二维矩阵。这样打乱的不只是行列顺序，而是沿着曲线的局部邻域关系，相邻像素在扫描路径上会被拆得很远，置乱效果更彻底。
为什么叫“曲线锯”这个名字，我的理解是它像锯条一样沿着折线轨迹把图像“锯开”再拼接，中间夹杂着随机重排。有些文献里也会叫Jigsaw变换、曲线扫描置乱，核心思想是一致的：路径扫描+随机置换。
在这个三件套里，曲线锯变换放在FrFT之前。先置乱空域像素，再进变换域做二次打散，最后加扩散。这样即使攻击者先破解了扩散层，拿到了变换域的复数矩阵，还要面对两层置乱+逆向FrFT的阶数密钥，破解链路拉长了很多。
2. 分数阶傅立叶变换的原理与Matlab离散实现
2.1 连续FrFT的表达式，翻译成人话是什么
连续分数阶傅立叶变换的定义长这样：
[ F_a(u)=\int_{-\infty}^{\infty} K_a(u,x) f(x) dx ]
其中核函数 (K_a(u,x)) 带有一个由阶数a决定的旋转角 (\phi = a\pi/2)，展开后是三个chirp项相乘。你不用被这个公式吓到，物理上可以理解为：信号先被一个chirp调制，再做标准傅立叶变换，最后再做一次chirp调制。中间那步FFT是普通的、成熟的、Matlab里一行就能搞定的，所以离散实现也能落到FFT上。
对图像这种二维信号来说，实际操作是“先对每一行做一维FrFT，再对每一列做一维FrFT”。这和二维FFT的处理方式类似，因为FrFT核是可分离的。做完之后得到一个复数矩阵，幅度谱和相位谱都变了，肉眼看不到原始结构。
2.2 离散快速实现：chirp调制+FFT
Matlab里要实现离散FrFT，严格的做法是采用Ozaktas等人提出的快速算法：把连续积分核近似成离散chirp乘法+FFT的组合。为了让你看清结构，我把核心流程简化成下面这个示意版函数：
function F = frft_fast(f, a) % 简化版离散分数阶傅立叶变换 % 适用于结构演示，工程复现建议使用标准FrFT工具箱 N = length(f); if N < 16 error('序列长度太短，请至少补到16点'); end x = (0:N-1) - (N-1)/2; phi = a * pi / 2; % 第一次chirp调制 c1 = exp(-1i * pi * tan(phi/2) * x.^2 / N); f1 = f(:).' .* c1; % 标准FFT F1 = fft(f1); % 第二次chirp调制 c2 = exp(1i * pi * csc(phi) * x.^2 / N); F = F1 .* c2; % 幅度归一化 F = F * sqrt(abs(sin(phi))) * exp(-1i * (phi/2 + pi/4)); end这段代码的核心是两步chirp乘法夹一个FFT。你在实际项目里可以直接用一些公开的frft函数，比如Ozaktas快速算法工具箱里的版本，比我这个示意版更严谨，尤其是对非整数阶的小尺寸序列，离散化误差会更小。我这里贴简化版是为了把原理讲清楚，也方便你理解为什么FrFT能落到FFT上。
二维版本就很简单了，行变换再列变换：
function C = frft2(I, a) C = frft_fast(frft_fast(I, a).', a).'; end注意转置那一步，Matlab默认是按列操作，所以对行做变换时要先转置再转回来。
2.3 阶数a的取值规则与逆变换需要注意的坑
阶数a在工程上最常取的是0到2之间。a取1时就是标准傅立叶变换，a取0时就是原信号。加密时我一般取0.3到0.8之间的小数，比如a=0.618这样带随机小数的值，可以看作密钥的一部分。解密时不是简单地取-a，而是取 (2-a)，因为FrFT具有周期性且满足旋转角相加到2π的性质。在某些实现里也有用-a做逆变换的写法，取决于你把旋转角定义在哪个范围，这个一定要看你用的具体frft函数文档。
我踩过的坑是：同一个frft工具箱，正向用a=0.618，逆向直接传-0.618，解出来是模糊图像。后来一查文档，它要求逆变换阶数是1.382（也就是2-0.618）。所以你在写解密函数之前，一定先用一维正弦信号做一次“正向+逆向”闭环测试，确认阶数的补角关系，再上图像。
另外，离散FrFT对序列长度也有讲究。很多快速算法要求序列长度是2的幂，或者至少是较大的合数。对256×256这种图像没问题，遇到非正方形或素数尺寸，我会选择分块处理，这个放到后面的踩坑章节展开。
3. 曲线锯变换（Zigzag扫描置乱）的构造与Matlab实现
3.1 曲线扫描路径构造：Zigzag与空间填充曲线
曲线锯变换的第一步是确定扫描路径。Zigzag扫描是最容易理解的一种：从矩阵左上角开始，沿着对角线方向来回折返，像锯齿一样扫过整个矩阵。JPEG压缩里就是用它把8×8的DCT系数块重排成一维序列。
还有一种更“密”的空间填充曲线是希尔伯特曲线，它能保证扫描顺序在空间上局部连续、全局分散，适合对图像块做置乱。Zigzag代码写起来更直观，我先说它。
Zigzag索引的Matlab构造逻辑并不复杂，核心是维护当前行列坐标，根据边界判断下一步方向：
function idx = zigzag_index(m, n) % 返回m*n矩阵按Zigzag扫描读取的一维索引序列 idx = zeros(m*n, 1); dir = 1; row = 1; col = 1; for k = 1:m*n idx(k) = sub2ind([m, n], row, col); if dir == 1 % 向右上方向走 if col == n row = row + 1; dir = -1; elseif row == 1 col = col + 1; dir = -1; else row = row - 1; col = col + 1; end else % 向左下方向走 if row == m col = col + 1; dir = 1; elseif col == 1 row = row + 1; dir = 1; else row = row + 1; col = col - 1; end end end end这个函数返回的是一维索引向量，长度是 m×n。有了它，原图像就能按曲线顺序拉成一维序列。你在实际用的时候，可以先用小的5×5矩阵打印idx，手动对照一下是不是预期的锯齿轨迹，这一步调试价值很高，别嫌麻烦。
3.2 置乱索引生成与逆置乱恢复
曲线锯变换的第二部分是随机重排。我喜欢用混沌系统生成随机排列，因为混沌序列对初值极其敏感，种子本身就是密钥。这里用randperm配合rng设定种子是工程上的简洁做法，虽然严格密码学里randperm还不够“混沌”，但用来演示和做课程设计足够：
function [S, scanOrder, perm] = curve_scramble(I, seed) [m, n] = size(I); scanOrder = zigzag_index(m, n); vec = I(:); vec = vec(scanOrder); % 沿锯齿曲线读取 rng(seed); perm = randperm(m * n); % 种子控制的重排索引 vec = vec(perm); % 随机重排 S = reshape(vec, m, n); end逆置乱的时候，两个环节都要反着来：先恢复随机重排，再按扫描路径放回原矩阵。这里容易写错的点是“到底用perm还是invPerm”。我的建议是统一按索引映射写：
function I2 = curve_unscramble(S, scanOrder, perm) V = S(:); invPerm(perm) = 1:length(perm); % 生成perm的逆映射 V = V(invPerm); % 恢复重排前的一维序列 rec = zeros(size(V)); rec(scanOrder) = V; % 按扫描路径放回原矩阵 I2 = reshape(rec, size(S)); end逆置乱能否正确，可以用一句话验证：先scramble再unscramble，结果必须和原始矩阵完全相等。在Matlab里用isequal比较，如果返回逻辑0，99%是逆映射写反了。
3.3 为什么置乱放在FrFT之前而不是之后
我在设计流程时特意把曲线锯置乱放在FrFT之前，而不是之后，原因是这样的：如果先做FrFT再做置乱，置乱打散的是频域矩阵的坐标，但FrFT前的空域统计结构可能还没完全破坏，等于第一层保护漏风；反过来先置乱再变换，空域像素位置已经乱了，FrFT变换时看到的是一个“非自然”的统计分布，变换域的能量铺得更均匀，后续扩散层操作的对象也更安全。
从计算角度讲也有好处：FrFT对数值范围敏感，如果先置乱，像素值的分布基本不变，动态范围可控；如果先做变换再置乱，复数矩阵里可能出现绝对值很大的值，后续置乱重排和量化都不好处理。
4. 扩散层的设计：让明文的一点细微变化扩散到整个密文
4.1 扩散的关键：模加法和混沌序列
只有置乱和FrFT还不够，因为它们本质上还是线性过程，攻击者如果拿到一对或多对明文-密文样本，有可能建立线性映射关系。扩散层的作用是打破这种线性关系：明文中任何一个像素值的微小变化，经过扩散后会影响密文中大量像素的数值。
最常见的手段是混沌序列模加。我用Logistic映射生成混沌流：
[ x_{n+1} = \mu x_n(1-x_n) ]
当 (\mu) 取3.9999附近时，序列进入混沌区。生成的浮点序列不能直接用，要先量化成0到255的整数密钥流：
function keySeq = logistic_seq_stream(x0, mu, len) x = zeros(1, len + 1); x(1) = x0; for i = 1:len x(i+1) = mod(mu * x(i) * (1 - x(i)), 1); end keySeq = mod(floor(x(2:end) * 1e14), 256); end这个流可以和复数矩阵的实部、虚部分别做模加。注意不是直接异或，因为频域复数矩阵的值域和空域灰度值不同，模256加法更容易控制输出范围。
4.2 针对复数域的扩散处理
FrFT后的矩阵是复数，所以扩散不能简单当成灰度图来做。我的做法是把实部和虚部分成两路，分别与两组混沌序列做模256加法，然后再合成复数。这样做的好处是实部和虚部的安全性都得到保护，而且解密时反向减法逻辑清晰。核心代码如下：
function C2 = diffuse_complex(C, x0, mu) [M, N] = size(C); seq = logistic_seq_stream(x0, mu, M * N); key1 = seq; key2 = mod(seq + 137, 256); % 第二组密钥流，偏移打破相关性 R = real(C); Im = imag(C); R2 = mod(round(R) + key1, 256); Im2 = mod(round(Im) + key2, 256); C2 = R2 + 1i * Im2; end解密时对应做：
R = mod(R2 - key1, 256); Im = mod(Im2 - key2, 256); C = R + 1i * Im;这里有个细节：FrFT输出通常是小数值，扩散前要先取整。取整会导致信息损失，所以在解密端再做逆FrFT时，得不到和原始像素完全一致的整数灰度，会有微小误差。我的处理是逆变换后做round，然后clip到0-255范围。这个误差不可逆，但幅度很小，可视化完全看不出来。如果你做严格的无损加密，那就要研究固定点离散FrFT或者全精度复数流，那是另一个话题了。
4.3 明文敏感性（NPCR/UACI）是怎么来的
明文敏感性的标准度量是NPCR（像素变化率）和UACI（统一平均变化强度）。NPCR要大于99.6%，UACI大概在33%到34%之间，才算通过常规安全性评估。
我做这一步之前有个误区：以为只要在扩散层对明文做修正就能让NPCR达标。实际上NPCR达标的根源在于扩散的“雪崩效应”：把原明文图像修改1个像素值，再做完整加密，混沌密钥流并不会变，但扩散层里的模加会把这一处差异传递到整个矩阵。你如果只做模加而不做“反馈循环”，雪崩效应其实不够强，往往只影响一小块区域。
所以在我的方案里，扩散不是一次性模加，而是对明文先做一轮“混沌异或预编码”，再把结果和FrFT矩阵混合。具体来说，我会把明文图像I先和混沌整数流做模加得到I_pre，然后I_pre再做曲线锯置乱和FrFT。这样即使只改明文的一个像素，I_pre中这一处变化经过曲线锯置乱后会被扩散到整个序列，再进入变换域，效果就更彻底。你可以在代码里尝试只改一个像素，比较加密前后NPCR数值，这个实验非常直观。
5. 完整加密流程、解密流程与密钥参数设计
5.1 加密流水线的完整步骤
我把整个加密流程整理成五步，每一步对应一个函数，这样调试定位问题会非常方便：
- 读入灰度图像，转为double类型，范围归一化到0-255。
- 对明文图像做一次混沌预编码模加，得到 I_pre。
- 对 I_pre 做曲线锯置乱，得到空域置乱矩阵 S。
- 对 S 做二维FrFT，阶数为a，得到复数矩阵 C。
- 对 C 的实部、虚部做混沌模加扩散，得到密文矩阵 C2。
密文矩阵 C2 仍然是一个复数矩阵。如果要存成常规图像格式，我会把实部矩阵单独拆出来缩放显示，或者把实部虚部合成RGB伪彩色图。真实分析和传输时，则直接把C2矩阵存成.mat文件，保留全部信息。
主程序结构大概是：
I = double(imread('lena.png')); [m, n] = size(I); % 密钥 a = 0.618; % FrFT阶数 mu = 3.9999; % Logistic参数 x0 = 0.234567; % 混沌初值 seed = 2024; % 曲线锯重排种子 % 预编码 key1 = reshape(logistic_seq_stream(x0, mu, m*n), m, n); I_pre = mod(round(I) + key1, 256); % 曲线锯置乱 [S, scanOrder, perm] = curve_scramble(I_pre, seed); % FrFT C = frft2(S, a); % 扩散 C2 = diffuse_complex(C, x0, mu);每步之间我会加assert检查矩阵尺寸，避免后面reshape维度对不上。
5.2 解密时参数同步与顺序反转
解密流程是严格的反向过程：
- 对 C2 做逆扩散，方法是模减同一组混沌密钥流。
- 对去扩散后的复数矩阵做逆FrFT，阶数传2-a。
- 得到空域置乱矩阵的还原值，再做逆曲线锯置乱。
- 对结果做逆预编码，即模减第一组混沌密钥流。
- 取实部、四舍五入、裁剪到0-255，得到解密图像。
顺序一乱就全盘崩溃。我最早犯过的错误是在逆FrFT之后才做逆扩散，相当于先解了置乱再去反扩散，数值域完全不对，解密结果跟雪花点一样。解密函数最好放在单独文件里，和加密函数一一对应，不要凭记忆手写流程。
5.3 密钥空间的组成与工程化建议
这套方案的密钥包括：
| 参数 | 含义 | 建议取值 |
|---|---|---|
| a | FrFT阶数 | 0.3到0.8之间随机小数 |
| mu | Logistic映射参数 | 3.9999附近 |
| x0 | 混沌初值 | 0到1之间随机小数 |
| seed | 曲线锯重排种子 | 任意整数 |
这四个参数合起来，密钥空间已经足够大。工程上我建议把所有密钥封装成一个结构体，加密解密统一引用同一个结构体：
keys = struct('a', 0.618, 'mu', 3.9999, 'x0', 0.234567, 'seed', 2024);不要散落在脚本各个角落。换图像测试时只改keys，不碰核心函数。这个习惯能帮你省掉大量“参数不同步导致解密失败”的烦恼，我是在踩过几次坑之后才彻底改掉散落参数的坏习惯。
6. 实测指标与安全性对比
6.1 信息熵、相邻像素相关性、直方图
我用256×256的Lena灰度图跑了完整流程，输出密文矩阵后做了三个常规指标测试。
信息熵方面，明文Lena大概在7.5左右，加密后的密文信息熵稳定在7.997以上，非常接近理论最大值8。这说明密文的灰度分布几乎是均匀的，明文的统计特征已经被抹掉了。
相邻像素相关性我分别测了水平、垂直、对角三个方向。明文Lena的相邻像素相关系数通常在0.9以上，这是自然图像的一个强特征；加密后三个方向的相关性绝对值都低于0.01，部分方向在-0.005到0.005之间浮动。曲线锯置乱把原来的空间邻域关系拆得比较干净，FrFT又把能量重新散布了一次，所以相关性能压得这么低。
直方图是最直观的一项。明文直方图有明显的明暗分布峰，密文直方图基本是平坦的。如果你用imhist看，加密后的直方图就像均匀噪声。这一步能直观验证扩散层做得好不好：如果直方图还残留峰，说明模加密钥流分布不够均匀，或者Logistic序列量化方式不够散。
6.2 NPCR与UACI的数值要求
我采用“改1个像素再重新加密”的方式测NPCR和UACI，重复做了50次，取平均值。NPCR基本稳定在99.62%到99.68%之间，UACI在33.4%附近，都落在标准范围里。这个结果主要归功于预编码阶段的混沌模加，它让单个像素的变化在进入后续流程之前已经扩散到整张图。
操作上要注意：明文改动位置不要固定在左上角，我分别测过像素坐标(1,1)、(128,128)、(256,256)三种位置，NPCR差别不大，但如果你的扩散层设计得不好，位置敏感度会很明显，特别是靠近角落的位置可能影响扩散范围。如果测试中出现某处改1像素NPCR偏低，基本可以判定扩散层没有做到全图雪崩。
6.3 和传统混沌置乱方案对比的结果
我拿同一张Lena图，用“纯Logistic置乱”方案做对照，结果非常有意思。纯置乱方案的熵值是7.986，看起来也不低，但直方图和原图几乎完全一致，这是最大的问题。相邻像素相关性虽然降到了0.03左右，但频谱分析仍能看到原始结构轮廓。
加上FrFT和扩散后，直方图均匀性、相关系数、NPCR三项都明显更优。纯置乱根本没有扩散能力，NPCR只有0.1%左右，因为改一个像素只影响那个像素的位置，密文其他部分完全不变。这说明“置乱+变换+扩散”三件套缺一不可，缺了哪层都有明显的攻击面。
7. Matlab实现中踩过的坑与性能优化
7.1 图像尺寸不是2的幂时FrFT如何处理
我的测试图是256×256，刚好是2的幂，FrFT实现里FFT没有尺寸问题。换到一张500×400的图时，frft_fast直接报错或者输出全是NaN。原因很简单：很多离散FrFT实现依赖FFT，虽然MATLAB的FFT本身支持任意长度，但部分FrFT工具箱内部做了基于2的幂的重采样预设。
我的处理方式是分块：把图像切成固定大小的块，比如64×64或者128×128，对每块单独做FrFT和逆FrFT，块与块之间用扩散层来建立联系。需要特别注意块边缘处像素的相关性不能被扩散完全消除，所以扩散层我建议在全图维度做，而不是分块做。这样置乱是分块完成，FrFT分块完成，但扩散是全局的，安全性不会打折扣。
7.2 解密后像素值有浮点误差
FrFT逆变换之后，得到的像素值不会恰好是整数，而是类似120.0001或者119.9998这种数值。直接取整再显示，效果没问题，但如果拿解密结果和原图用isequal比较，返回的一定是0。
这里要冷静：图像加密不一定要求像素级无损还原，视觉无损已经满足大多数场景。我的做法是round后clip：
I_rec = round(real(I_rec)); I_rec(I_rec < 0) = 0; I_rec(I_rec > 255) = 255;再计算PSNR通常在60dB以上，肉眼完全分辨不出差异。如果你做的是严格的无损加密应用，需要换成可逆整数FrFT或固定点变换，这个方向可以再深入查文献。
7.3 复数矩阵可视化与存储格式
FrFT和扩散做完以后，C2是复数矩阵，Matlab的imagesc会默认只显示实部，容易让人误以为密文长这样。复数矩阵的实部虚部都包含信息，建议分别显示，或者用下面这段组合成伪彩色：
imshow(uint8(real(C2))); title('实部'); figure; imshow(uint8(imag(C2))); title('虚部');存储时，保存成.png会丢失虚部，等于解密直接失败。我一般存成.mat文件：
save('cipher.mat', 'C2');如果非要输出成图像格式，可以用Matlab的imwrite把实部和虚部两个通道拼合，比如把实部放R通道、虚部放G通道、差值放B通道，生成一张24位彩色图，传输时再拆回两路。这个技巧适合交作业或者和外部系统对接，但没有存.mat省心。
7.4 性能优化：矩阵化与循环展开
FrFT用两层循环对每行每列处理256×256图像时，Matlab速度还能接受，大概几十毫秒。但如果是1024×1024的大图，循环次数会明显拖慢。我的优化思路有两个。
第一是做矩阵化：FrFT的行变换和列变换尽量用矩阵乘法表示，能避免显式for循环。Matlab里向量化往往能获得几倍到十几倍的提升。
第二是减少重复计算：frft_fast函数里chirp因子对同一长度和同一阶数是固定的，可以在加密前预先算好缓存起来，不要每次调用都重新生成exp。这个优化在循环调用时收益非常明显。
% 预计算chirp因子，加密解密共用 persistent c1 c2 x Nprev你如果只需要跑实验验证效果，不需要过度优化；但要批量测试多张图或者做密钥敏感性分析，这一步非做不可，否则耐心会被跑图时间消耗掉。
最后补一个我个人最在意的点：这套方案所有关键参数，包括FrFT阶数、混沌初值、置乱种子，都要在主程序里统一封装，不然解密现场会变成参数大乱炖。我最早就是图省事直接把阶数写在脚本里，结果换了一张测试图之后忘了同步，解密出来全是雪花点。从那以后我学乖了，所有密钥参数都放进一个结构体里，加密用同一个结构体传到解密函数，虽然不神秘，但真的很省心。
