比起隐写，我认为DFT更像是盲水印。
在各种D*T中，DFT看起来应该是相对简单的一个。DCT我有个大概印象，但DWT我实在不知道是什么高端的玩意儿。总之，本文将讨论在DFT变换域隐写数据。
什么是DFT
DFT全称Discrete Fourier transform，即离散傅里叶变换。那什么是傅里叶变换呢？我不知道该怎么给你解释，但你可以参考这个视频。简单地说，傅里叶变换是一种拆解信号的过程。以音频信号为例，一段录音对于计算机来说就是麦克风在时间尺度上对气压进行采样，一段录音可以看作是气压在时间上的函数。可是对于人的听觉来说，气压并不是特别有效的信息，而气压的变化才是。气压在某个固定周期内重复波动，这近似就算是频率了。我们的耳朵能够感受到频率，但计算机只能收录气压的采样，那怎么才能让计算机分析出人能够感知的频率呢？傅里叶变换就能够做到这件事。对于一维信号（只有一个变量），傅里叶变换将信号拆解成一系列正弦或余弦信号的叠加，而这里的正弦或余弦信号就是纯粹的频率。换句话说，傅里叶变换允许你用一系列纯粹的频率来表示一段音频（例如三倍的440Hz加上五倍的320Hz）。类似地，离散傅里叶变换是作用在离散信号上的。在理想情况下，声音应该是气压在时间上的连续函数，但对于计算机来说，声音则是每隔一段时间对气压的采样，它并不连续。因此将傅里叶变换推广到这种离散信号上，就变成了离散傅里叶变换。这在计算方式上变了，但在原理上没有改变。
对于图像这样的二维信号（两个变量，x、y），傅里叶变换同样将其拆解为频率，但这里的频率并不是人耳能听到的频率了，它指的是像素明暗的周期性变化，而这里的频率也具有两个参数，分别表示水平频率（x方向）和垂直频率（y方向）。如果一个频率是(5,0)，那么这个频率看起来应该是水平方向的明暗交替；而频率(0,5)则是竖直方向的明暗交替；而频率(5,5)则是从左上角到右下角方向的明暗交替。如果光看文字难以想象的话，我推荐观看这个视频。
从DFT到FFT
离散傅里叶变换的计算其实不难，并没有什么特别高神的东西：
$$ F(x)=\sum_{n=0}^{N-1}e^{-i\frac{2\pi}{N}nx}f(x) $$
其中F(x)是傅里叶变换，f(x)是原来的信号。如果是二维的话，根据二维傅里叶变换的分离性，可以分解为在两个方向上分别进行傅里叶变换：即先以行为单位，计算x方向的傅里叶变换，然后以列为单位，计算y方向的傅里叶变换：
suspend fun ComplexChannel.fft2d(): ComplexChannel = coroutineScope {
// result: [height, width]
val rowTransformed = this@fft2d.indices.map { y ->
// feed each horizontal line into fft
async { FFT.fft(this@fft2d[y]) }
}.awaitAll()
// here we calculate based on columns
// aka the output shape is [width, height]
val columnTransformed = this@fft2d[0].indices.map { x ->
// we need to collect each vertical line
async {
val columnData = Array(this@fft2d.size) { y -> rowTransformed[y][x] }
FFT.fft(columnData)
}
}.awaitAll()
// recover the structure back to [height, weight]
Array(this@fft2d.size) { y ->
Array(this@fft2d[0].size) { x ->
columnTransformed[x][y]
}
}
}但如果你仔细观察DCT的定义的话，你会发现这个算法其实并不高效。在处理$e^{-i}$那部分时，需要使用欧拉公式将其变换为三角函数，这对于计算机来说本身就是一个难事了，更加雪上加霜的是，计算傅里叶变换的一个点就要对整个输入空间遍历一次。我们先看按行计算的傅里叶变换，假设图片宽度是w，那么计算一个点的傅里叶变换就要遍历w次，而完整计算一行就需要遍历$w^2$次，这还是一行，如果图片的高为h，那么光是按行计算就要遍历$hw^2$次；接下来再考虑按列计算，类似地，要$wh^2$次遍历。当你看到复杂度里面带着平方的时候，大概就能想到它算起来不会很快。
为了解决这个问题，人们想出了FFT，快速傅里叶变换。由于我不是很在行数学，因此我就不多说其中的原理了。这里我参考了一篇文章，文章中提供了C++实现的能处理任意长度输入的FFT算法。在Java中，Apache commons Math3这个库提供了FFT实现，但它要求输入大小必须是2的幂。我认为padding不够优雅，且在内存上不够高效，因此我仿照文章中提供的C++实现编写了一个Kotlin下的FFT算法。
在变换域隐写的问题
目前为止，我们解决了计算的问题，接下来我们得想办法在频域中隐写。与LSB不同，我们之所以能够在LSB中隐写数据，主要有以下依据：
- 隐写的数据能够按照原样写入、读出，即写入是什么，读出还是什么
- 隐写的数据对原图像没有太大的破坏
而频域隐写则无法依赖这两条。傅里叶变换是浮点数运算，在修改频域后还需要使用对应的逆运算还原回图像，在这个过程中，由于浮点数精度受限，写入的数据并不保证会被原封不动的保存。此外，由于我们修改的是频域，因此任何一个点的修改都会影响到图片全局，尤其是人眼对低频信号更为敏感，因此隐写在低频部分会使得图像面目全非；可是如果隐写的频率太高，那么受限于采样定理（奈奎斯特定理），经过逆变换得到的图片可能难以保留我们写入的数据（例如图像的宽度为w，但隐写数据出现在超过w/2的频率以上，此时可能会出现采样定理所说的混叠和失真）。
除了以上问题，傅里叶变换的结果也不太利于隐写：通过将变换结果输出成CSV（输出复数的模），我观察到对于较大的图片（例如边长在1000像素以上的），在边缘部分（低频）的变换结果通常比较大，有时会达到两三万的程度，而中心部分（高频）则通常较弱，只有几百甚至几十，如此之大的差距使得变换结果没办法被线性地展示出来，同时，如此宽阔的取值范围也给隐写带来了很大的难度。
这里需要补充一点，我认为这一点隐含在了傅里叶变换中，但如果你没有那么深入理解这个变换的话，这一点补充或许能够解答你的疑惑。
在傅里叶变换的计算中我们同时使用了正弦和余弦，对于给定频率，都会得到两个分量幅度，这也就是为什么傅里叶变换的结果是个复数。如果以余弦信号来表示原始信号，那么这个复数的模就是余弦频率的强度；而复数的辅角（argument）就代表了这个频率的初始相位。
在图片处理中，频率可以看成是一个纹理，和一维类似，强度决定了这个纹理在图片中的占比，而相位则决定了这个纹理在图片中的初始位置。如果保留强度而调整相位（调整实部和虚部的比例），则相当于改变一个纹理的位置；而保持相位改变强度（等比例同时缩放实部和虚部）则相当于调整这个纹理的明暗。本文只修改纹理的明暗，不改变纹理的位置。
此外，由于傅里叶变换的对称性，在修改时需要保持共轭，这一点会在后续展示代码的时候强调。至于为什么傅里叶变换会有对称性，为何这种对称性要求复数共轭，可以参考这篇知乎回答。
对数可视化与盲水印隐写
变换结果的跨度太大，那我们就把跨度缩小好了。对数正好可以用来缩小跨度：
$$ L(v) = \log_{10}(1 + k \cdot v) $$
式中，v是复数的模，k是一个可变的系数，用来调整增益，这样一来函数L的输出就变成了正实数，如果再将其输出裁剪到0到1之间，那就可以产生一副灰度图了。
既然都有图了，不如我们把思路打开：如果频域隐写没办法忠实地还原我们写入的数据，那让数据能够容错就是了。此外，对于图片的旋转、裁剪等也会反映到频域上，因此我们需要一种能够在旋转、裁剪、缩放的情况下还能解码的编码方式。听起来是不是有点耳熟？这不就是二维码吗？
DFT隐写常常作为盲水印使用，因为它的编码效率不高，但是能保证你隐写进去的东西足够健壮，不会因为平移、旋转、裁剪、遮挡等攻击而消失不见。如果我们将数据编码成二维码，然后将二维码作为图片数据嵌入到频域，这样无论具有怎样的失真、旋转、拉伸，只要还在二维码的容错范围内，我们就可以读出其中的数据。
不过这种隐写方式有一个弊端：那就是不够隐蔽。在LSB隐写中，我们可以用加密、随机等方式将数据分布到不同的通道，这是因为我们的编码方式依赖绝对位置。而DFT的一个特征就是频域能够随着图片的旋转而旋转，即便图像经过剪裁，也能在一定程度上保留频域数据，这就意味着DFT隐写时，我们不能依赖绝对位置，否则我们隐写的数据就和LSB一样脆弱了。可是没办法依赖绝对位置，也就没办法用随机性隐藏数据，尽管我们可以对数据加密，但频率里一个工工整整的二维码，这多少有些可疑的离谱了。但是没办法，你获得了健壮性，就要牺牲隐蔽性。如果你要隐蔽性，那就只能牺牲健壮性。
关于这种隐写的容量，那肯定是不如LSB，但经过测试，用这种方法传递ed25519的公钥还是没问题的。如果你的图片选用得当，我认为基本的文字传递还是可以做到的。
黑白半幅傅里叶隐写
综上所述，我们可以将黑白图片编码到频域中，而取决于我们使用的象限，我们可以：
- 使用一二象限，将数据编码在频域的上半部分，而下半部分是对应的共轭
- 使用二三象限，将数据编码在频域的左半部分，而右半部分是对应的共轭
在开始写代码之前，老规矩，先写接口：
interface Fourier {
suspend fun encode(
source: BufferedImage,
redMask: BufferedImage? = null,
greenMask: BufferedImage? = null,
blueMask: BufferedImage? = null,
): BufferedImage
suspend fun decode(
input: BufferedImage,
energyLogKList: List<Double>
): List<Triple<BufferedImage, BufferedImage, BufferedImage>>
suspend fun decode(
input: BufferedImage,
energyLogK: Double
): Triple<BufferedImage, BufferedImage, BufferedImage> = decode(input, listOf(energyLogK))[0]
}
这里规定了编解码的函数签名。注意到解码的时候需要提供一个energyLogKList，这个就是刚才提到的对数变换的K的取值列表。一般编码的时候我们是根据给定的K来计算的，但在可视化的时候可以一口气产生不同K值下的可视化结果。虽然我们可以直接返回傅里叶变换的结果，但我认为这部分属于是内部实现，我不希望使用者获知其中的具体细节。你可以观察到，无论是编码还是解码，其中设计的数据类型全都是Java标准库自带的，并没有涉及到其他的库（例如Apache commons Math3的Complex类）。
前述两种编码，无论使用哪种方式，编码的区域是不变的。因此我们只需要让子类给我们画出一片水印覆盖的区域就可以在父类中做到对称修改频域了。
abstract class BlackAndWhiteHalfFourier(
private val targetAmp: Double
) : Fourier {
protected abstract fun checkInputSize(width: Int, height: Int, mask: BufferedImage)
protected abstract fun calculateBox(
width: Int, height: Int,
mask: BufferedImage
): Pair<Pair<Int, Int>, Pair<Int, Int>>
private fun encodeChannel(channel: ComplexChannel, mask: BufferedImage) {
val height = channel.size
val width = channel[0].size
checkInputSize(width, height, mask)
val (p0, p1) = calculateBox(width, height, mask)
val (x0, y0) = p0
val (x1, y1) = p1
// apply encoding
for (y in 0 until mask.height) {
for (x in 0 until mask.width) {
val oldComplex = channel[x0 + x, y0 + y]
val targetAbs = 9.0 / targetAmp
val oldAbs = oldComplex.abs()
val newComplex = when (mask.getRGB(x, y)) {
Color.BLACK.rgb -> Complex.ZERO
Color.WHITE.rgb -> if (oldAbs < targetAbs) oldComplex * (targetAbs / oldAbs) else oldComplex
else -> continue
}
channel[x0 + x, y0 + y] = newComplex
channel[x1 - x, y1 - y] = newComplex.conjugate()
}
}
}
private fun encodeChannelOrNull(channel: ComplexChannel, mask: BufferedImage?): ComplexChannel =
channel.also { if (mask != null) encodeChannel(it, mask) }
final override suspend fun encode(
source: BufferedImage,
redMask: BufferedImage?,
greenMask: BufferedImage?,
blueMask: BufferedImage?
): BufferedImage = coroutineScope {
val fft = source.fft2d()
val redChannel = async { encodeChannelOrNull(fft.first, redMask) }
val greenChannel = async { encodeChannelOrNull(fft.second, greenMask) }
val blueChannel = async { encodeChannelOrNull(fft.third, blueMask) }
ComplexImage(
redChannel.await(), greenChannel.await(), blueChannel.await()
).ifft2d().recoverToImage()
}
final override suspend fun decode(
input: BufferedImage,
energyLogKList: List<Double>
): List<Triple<BufferedImage, BufferedImage, BufferedImage>> = coroutineScope {
val recoveredFFT = input.fft2d()
energyLogKList.map { k ->
val red = async { recoveredFFT.first.visualize { it.energyLog(k) } }
val green = async { recoveredFFT.second.visualize { it.energyLog(k) } }
val blue = async { recoveredFFT.third.visualize { it.energyLog(k) } }
Triple(red.await(), green.await(), blue.await())
}
}
}这里省略了不少工具类。可以移步Github仓库查阅完整代码。
这里我们按照RGB三个通道分别编码。为了加速计算，这里使用了Kotlin协程。编码时取决于对应通道的mask是否为null，如果为null则不编码，否则将mask编码到对应通道中。编码时让子类给出一个中心对称的区域，父类中只要遍历mask，对称地修改这个区域内的频域即可。解码时直接将对应的通道做傅里叶变换，然后可视化即可。
对于子类，以使用一二象限的编码方式为例：
class BlackAndWhiteHorizontalFourier(targetAmp: Double) : BlackAndWhiteHalfFourier(targetAmp) {
override fun checkInputSize(width: Int, height: Int, mask: BufferedImage) {
require(mask.height < height / 2) { "Mask too big" }
}
override fun calculateBox(width: Int, height: Int, mask: BufferedImage): Pair<Pair<Int, Int>, Pair<Int, Int>> {
// top left
val x0 = width / 2 - mask.width / 2
val y0 = height / 2 - mask.height
// bottom right
val x1 = width / 2 + mask.width / 2
val y1 = height / 2 + mask.height
return (x0 to y0) to (x1 to y1)
}
}在某种程度上来说，calculateBox的计算结果也是共轭的。类似地，这是使用二三象限的编码方式：
class BlackAndWhiteVerticalFourier(targetAmp: Double) : BlackAndWhiteHalfFourier(targetAmp) {
override fun checkInputSize(width: Int, height: Int, mask: BufferedImage) {
require(mask.width < width / 2) { "Mask too big" }
}
override fun calculateBox(width: Int, height: Int, mask: BufferedImage): Pair<Pair<Int, Int>, Pair<Int, Int>> {
// top left
val x0 = width / 2 - mask.width
val y0 = height / 2 - mask.height / 2
// bottom right
val x1 = width / 2 + mask.width
val y1 = height / 2 + mask.height / 2
return (x0 to y0) to (x1 to y1)
}
}编码结果
以下展示的全都是截图，因为原图太大了。
这是编码之后的图片，肉眼察觉不出异常。
这是解码得到的：
可以看到这张图片并没有太多细节，所以中心部分的高频比较弱。但我们的Data Matrix中包含了白底来增强对比度，因此使用对应的app可以无压力的扫描。接下来我们试试看用于分析LSB隐写的方法。
直方图上看不出明显的异常，因此这里就不放了，这里主要看位平面分析和随机色彩映射：
左边是隐写后的图片，右边是原图。可以看到原图的白底部分在位平面上非常干净，而经过隐写后的位平面则充满噪点。随机色彩映射也能看到类似的问题。但和LSB隐写引入的噪点明显不同：
右侧是LSB隐写的分析结果，可以看到DFT隐写带来的噪点看起来更柔和一些，例如树叶边缘看起来就像是在发光一样，而LSB隐写带来的噪音则相对生硬。当然，正如我们在上一篇文章所说，噪点不是问题，只要我们的图片有充分存在噪点的理由，那就没有异常。
接下来我们看看DFT隐写的能力。前面说了半天什么拉伸啊旋转啊，说了一大堆，我们不如亲眼看看。当然了，看二维码旋转没有意思，我准备了一大段文字作为mask来隐写。这里依次对隐写结果做缩放（等比例和非等比例）、旋转和遮挡，值得一提的是这些处理都保持了原有的尺寸，这是因为缩小尺寸将导致DFT无法计算高频部分。
首先是原图：
上半部分是肉眼可见的图片，下半部分是DFT后得到的结果。可以看到文字能够相对清晰的被识别，而上方的图片也没有过于明显的变化。
接下来我们看看等比缩放和平面位移：
根据傅里叶变换的性质，无论图片如何位移，傅里叶变换都能提取出其中的频率信息，但是这里进行了等比缩放，可以看出，水平方向上因为隐写数据躲开了高频的位置，经过缩放后虽然叠在一起，但通过对称性能够从右侧补全数据。而垂直方向就没有那么幸运了：垂直方向的数据一部分分布在中心的高频位置上，经过压缩后高频信息被破坏，因此文字的结尾几行已经无法辨认了。
接下来我们看看非等比缩放加上10度的旋转：
因为旋转的原因，DFT的结果也跟着旋转了10度，由于这次缩放主要发生在水平方向，因此垂直方向的数据还能够被辨认。通过拼凑左侧和上侧的图案，除了最后一行之外，大部分文字还是可以阅读的。
接下来看看部分数据丢失会怎样。这里对原图的四边进行了裁剪，并拉伸放大回原来的分辨率：
看起来部分数据的丢失并不影响隐写，这是因为隐写发生在频域，而频域对图片的影响是均匀发生在所有地方的。可以看到隐写的文字依旧清晰可见。
当然，光是裁剪还不够，我们直接在图片上涂鸦：
如果是LSB编码，这种涂鸦基本上就会破坏掉相当一部分数据了。如果把涂鸦换成白噪声（这里没有演示），那对于LSB来说就是灭顶之灾了。但是这种简单的涂鸦在频域上属于低频，因此对隐写数据的影响有限。如果换成噪声的话，取决于噪声的频率，可能会影响隐写数据。
最后我们来看看两个更加过分的操作。第一个是假如半径为2像素的高斯模糊：
可以发现图片看起来稍有模糊，但高斯模糊平滑了像素，拉低了频率，使得我们的隐写数据完全被破坏了。
最后我们看看JPEG压缩：
这里使用了两种不同的压缩质量。上面是压缩质量8，在Photoshop里算是中等质量，下面是压缩质量12，在Photoshop里算是最佳。可以看到经过JPEG压缩后，虽然能够隐约看出隐写数据，但高频部分几乎完全不可读了。
这里使用的mask其实有点不完美，因为图片比较宽，导致隐写的数据在水平方向上的频率比较低，而垂直方向则集中在相对高频的部分，这部分数据受到挤压、裁剪等动作后也会随之丢失。在实际编码中，应该多加测试，保证隐写的数据躲开易受破坏的高频，同时保证尽可能不影响图像（主要是低频）。
结语
本文尝试了DFT隐写。正如开篇所说，我认为DFT在隐写方面并不隐蔽，一个简单的傅里叶变换就足以让肉眼观察到可疑的隐写迹象。也正因为这个原因，DFT更多被用在盲水印追踪领域，这个领域和隐写不同：隐写旨在藏匿数据，首要指标是隐蔽性和较大的数据量；而盲水印则通常只是隐写少量能够表明版权归属的文字或图案，同时追求健壮性，即能够抵御各种图像处理手段，经过处理后仍然能解读出水印以供查证版权。
这个系列的下一篇文章要是不出意外的话，将会写DCT（离散余弦变换）隐写。DCT在JPEG压缩中广泛使用，但我并不打算做成JPEG隐写，因为那主要是针对文件格式进行的（修改系数或量化表）。在搜索相关文献的过程中，我看到一种DCT-SVD的隐写方法，该方法对DCT系数进行奇异值分解（SVD），然后奇异值中藏匿数据。虽然我还不知道该怎么做，但是我十分好奇这个过程，也打算尝试一番。
-全文完-
【代码札记】初探图片隐写术（Steganography）之DFT隐写 由 天空 Blond 采用 知识共享 署名 - 非商业性使用 - 相同方式共享 4.0 国际 许可协议进行许可。
本许可协议授权之外的使用权限可以从 https://skyblond.info/about.html 处获得。