信号图像处理：傅立叶变换与频域分析
最近在QQ群里看到了下列一则消息：
豆瓣在网页和APP中加入盲水印 可追踪截图用户
豆瓣在15日于其网页版标题下方加入了一行小字，内容是uid、tid和截图时间。这些文本颜色与网页背景色大体接近，肉眼无法直接分辨，如果有人截图通过调节对比度就可以显现出来。豆瓣可以轻松追踪截图者是谁。
在被曝光后豆瓣已将该部分设置为禁止选中。另外豆瓣早在2021年已经在其App中加入了类似功能。
隐写技术并不是新鲜事，其它互联网公司也早有应用，甚至更加先进。
在此前阿里巴巴月饼事件中泄露内网截图的员工虽然对图片做了处理，但是依然可以通过频域手段添加的数字盲水印精确定位到泄密员工。知乎用户 fuqiang liu 经过测试，使用涂抹，剪切，放缩，旋转，压缩，加噪，滤波等手段无法过滤这种盲水印，唯一的办法是对屏幕拍照。但是，依然可以通过文本间距等手段区分不同用户。审查者在这方面有很多技术选项。
有不少同学十分好奇其中提到的“通过频域手段添加的数字盲水印”是如何实现的。对于初学者而言也确实很难将一张图片与所谓的频率联系起来，那么到底什么是数字图像的频域分析呢？
傅立叶变换
引用一段Wiki对于傅立叶变换的介绍：
傅里叶变换（法语：Transformation de Fourier、英语：Fourier transform）是一种线性积分变换，用于信号在时域（或空域）和频域之间的变换，在物理学和工程学中有许多应用。因其基本思想首先由法国学者约瑟夫·傅里叶系统地提出，所以以其名字来命名以示纪念。实际上傅里叶变换就像化学分析，确定物质的基本成分；信号来自自然界，也可对其进行分析，确定其基本成分。
可知，傅立叶变换做的工作就是把一段信号从时空域搬到频域上，频域自然就是指根据不同频率分布的数值，那么这个频率到底指什么的频率？傅立叶级数告诉了我们如下的答案：
在对傅里叶级数的研究中，复杂的周期函数可以用一系列简单的正弦、余弦波之和表示。傅里叶变换是对傅里叶级数的扩展，由它表示的函数的周期趋近于无穷。
wiki上用了一张图相当形象地展示了这一过程：
即，我们将以时间变换的x轴替换为了根据频率变换的x轴，而同一个三角函数的频率是不变的，所以频域空间上得到的是一条条线段。
傅立叶也为我们提供了一个方便的公式来进行傅立叶变换（连续），对于可积分的函数\(f\)：
\[\hat{f}(\xi)=\int^\infty_{-\infty}f(x)\exp{(-2\pi ix\xi)}\text{d}x\]
这样，我们可以再通过欧拉公式把其变为正弦、余弦函数的组合，即做到了将所有的函数都转换到频域上的需求。
离散傅立叶变换
可是如果仅仅知道连续傅立叶变换，并不能很好地解决我们的问题，毕竟图像由一个个像素组成，属于离散信号，这个时候应该怎么做傅立叶变换呢？很简单，类似于黎曼积分，对于\(N\)点序列\(\{x[n]\}_{0\le n \le N}\)，离散傅立叶变化为：
\[\hat{x}[k] = \sum^{N-1}_{n = 0}\exp{(-i\frac{2\pi}{N}nk)x[n]}\]
图像的傅立叶变换
虽然傅立叶变换看起来很美好，可以描述任意信号的频域分布，但是一张图片不能是一条线组成的啊，这里就需要用到二维的傅立叶变换了，所谓二维傅立叶变换，也不过是在空域的\(x\)轴再加上一个\(y\)的维度，这里我们只讨论离散情况：
\[F(u, v) = \sum^{M-1}_{x=0}\sum^{N-1}_{y=0}f(x,y)\exp{(-j2\pi(ux/M+vy/N))}\]
\[\exp{(-j2\pi(ux/M+vy/N))} = \cos\theta + j\sin\theta, \theta = 2\pi(ux/M+vy/N)\]
那么应该怎么确定一个正弦/余弦平面波？首先肯定的，有三要素：频率\(\omega\)，幅度\(A\)，相位\(\phi\)，同时，平面有着一维函数不具有的性质：方向\(\vec{n}\)，即两个平面波可以拥有相同的三要素，但是却像海浪一样，拥有不同的波动方向。
K空间
将一个二维函数分解为平面波之后，我们应该怎么形象地可视化这些波的不同分布呢？已知每个波有不同的法向量\(\vec{n}\)，频率\(\omega\)，那不如就用这两者构建一个空间，这个空间中的一个坐标点\((u,v)\)代表了在法向量\(\vec{n} = normalize(u, v)\)方向下，频率为\(\sqrt{u^2 + v^2}\)的一个波，每个波用0-1表示幅度，即可构建一个k空间矩阵，用来存储各个波的分布，这样，我们就可以很自然地得到一个中心代表低频、四周代表高频、不同方向代表不同法向量、不同明暗代表不同幅度的图像（这里我们选择性忽略了相位）。顺带提一句，磁共振也就是利用这个原理，先构建k空间中的信息，再利用逆向傅立叶变换，得到原本的医学影像的。
离散二维傅立叶变换的代码实现
根据上面的公式，我们就可以编写出如下的离散二维傅立叶变换的代码了，这里为了方便，选用了Python作为我们的编写语言，当然相应的，执行速度非常的慢，因此选择将图片resize到原来的0.05倍大小来查看：
1 | def process_img(raw_img): |
可以看到很简单地，几乎和公式保持完全一致地求和，这里的对频域图片的每个像素求值的时间复杂度达到了恐怖的\(O(n^2)\)，但是可以看到，效果还是非常不错的（原图片-自己编写的dft2-np提供的fft2）：
当然，这里最后的图像是将(u = 0, v = 0)移动到中心后得到的图像（np.fft.fftshift）。从得到的dft图像中我们得出这张肥猫小猫的图片的几个特征： - 在(0, 0)的地方振幅几乎为0，周围亮度高：超低频缺失，几乎不存在一成不变的区域，反之我们可以构造一副纯白的图片，其在(0, 0)点处应该为白色（下图需要放大才能看清那个白点）
- 频域图像有黑色噪声分布：图像本身有不连续的噪声出现，是缩放导致的图片噪声，将图片放大后消失
图像频域分析的应用
说了这么多，频域上的分析到底有什么意义呢？平时我们很多直接分析图片发现不了的信息，都可以统统甩到频域上解决，或者是从频域的角度出发，对图片进行修改，比如我们已经知道了，高频信息往往出现在图片中的边缘区域（临近区域变化明显），这个时候，便可以通过过滤掉这些频率的信息（高斯模糊滤波）来实现；相反的，如果我们需要做边缘检测，便可以提取出频域上的高频区域来进行边缘的检测。
如何在一张图片中隐藏信息
在这方面其实有很多研究了，比较简单的有LBS（最低有效位），通过使用图片的每个像素最低bit来实现信息的嵌入，但是很明显，这种方法仅仅适用于无损图片，图片经过压缩或者变换之后很容易导致信息的缺失，这对于读准确率需求不高的情况是一个优秀的解决方案；对于需要高读准确率的场景，一方面会造成信息的浪费，另一方面在较高信息率的情况下，纠错正确率很难得到保证，这个时候，频域隐写技术就是一个很好的替代品。
频域变换技术通常使用的是不同于上文的DCT（离散余弦变换），DCT 变换的好处是，如果原序列是实数序列，那么变换后也是实数序列。[1]这里我们简单了解一下DCT隐写算法。
DCT隐写（两点法）
DCT隐写算法直接针对JPEG图像进行（JPEG图像天然以8x8的子块存储），隐写算法的过程如下：
- 取出一个8x8的分块
- 通过DCT将分块变换到频域空间
- 随机（全局一致）选取两个点，并且根据需要嵌入的比特，通过像素值值赋予这两个点顺序（如果嵌入0则这两个点必须小的在前面，否则大的在前面，通过交换这两个像素来实现）
- 通过IDCT将块变换回去
这里需要注意的是，如果希望信息的抗干扰性高一些，我们需要在选取两个点的时候使其距离相对比较大。
提取的算法就很显而易见了：
- 对于每个8x8的分块做DCT
- 根据值大小，检查是否存在交换的现象，提取该分块的比特值
- 合并信息
总结
这一篇博客算是我写给自己的一个交代，之前一直无法理解一张图片是如何和信号、频域关联起来的，最近几天花了不少时间，查阅了大量的资料，总算是有了一个不算模糊的认知。
频域的分析也是实时和离线渲染中常常使用的工具，特别是实时渲染的后处理，经常通过频域搞事情，故希望自己在遇到这些东西之前，把图像频域方面的事情搞明白。
参考
信号图像处理：傅立叶变换与频域分析