喜欢科普的我思考良久，准备挑战一下这个高难度的问题。难并不是因为有多么难学，而是非常的抽象，无法直观理解，也是很多在校大学生的梦魇。借用我在《线性代数有什么用？学习线性代数的意义在哪？》里写过的一句话：将一门知识讲复杂了很容易，但是讲容易了的确是非常难的，这需要极高的学识和深刻的理解。相信放在“傅里叶变换”中也是再合适不过的了。下面进入正文。傅里叶变换的核心是从时域到频域的变换，而这种变换是通过一组特殊的正交基来实现的。这里面涉及到很多专业术语，首先来解释为什么一定要进行时域到频域的变换呢？时域里有哪些我们解决不了的问题吗？时域时域是描述一个数学函数或物理信号对时间的关系，这也是我们日常中最容易直观感受的一种域。从我们学物理开始，很多物理量的定义都是跟时间相关的。速度：位移与发生这个位移所用的时间之比电流：单位时间里通过导体任一横截面的电量功率：物体在单位时间内所做的功的多少...很多物理量的定义都是基于单位时间产生的效果或者变化，以时间为参考让我们更容易理解。但是容易理解不代表方便使用，或者说方便计算。比如我截了一段音频的波形图（来自李荣浩《麻雀》的副歌部分——“我飞翔在乌云之中，你看着我无动于衷...”），如下图。波形图其中横轴是时间t，纵轴是振幅A [-1, 1]。假设播放器读入这段音频进行音频播放。现在我想让音量大一些，播放器应该怎么做呢？因为上面的波形图的振幅对应的其实就是声音的强度，如果想让音量大一些，只需要将整体的振幅同比例扩大即可。这个需求看起来很容易满足。但如果我比较喜欢低音效果，虽然李荣浩的音色已经比较低了，但是我还是想加强上面这段音乐的低音部分，更加厚重一些，那播放器应该怎么做呢？虽然这是一段美妙的音乐，但是从时域的图像看起来，似乎杂乱无章，想找到低音部分根本无从下手，跟不用说将低音部分加强了。因为高中低音在时域中是杂糅在一起的，我们无法将他们剥离开来，随便改动波形图中的一小部分，都会同时影响到高中低音。所以如果播放器仅仅对时域信号进行处理是无法完成这个需求的。和时域的这种限制类似的还有RGB空间。任何一个颜色都可以通过R/G/B(红/绿/蓝)三原色表示出来。如下图。三原色通过调整三种颜色的配比，就能混合出各种颜色。为什么我们经常通过RGB空间来表示所有的颜色呢？因为人类有三种视锥细胞，而这三种视锥细胞最敏感的波长接近于红/绿/蓝（如下图）。所以任何颜色对于大脑来说，都是这三种视锥细胞电信号的混合作用。这也是我们使用RGB空间的生物学基础。三种视锥细胞虽然RGB空间和我们的视锥细胞原理类似，而且模型非常简单。但是在某些条件下，它仍然无法满足我们的需求。比如，我们在拍照时有时会出现红眼现象（如下图的美女）。红眼我们需要PS掉红眼，但是我们如何在RGB空间中找到红色的范围呢？有人可能会说，R值越大的地方代表越红，是这样的吗？我们看(R,G,B)=(170, 0, 0)时，颜色如下，#AA0000上图的颜色我们可以认为是红色的范围。但当RGB=(187, 187, 187)的时候，颜色如下图所示，#BBBBBB虽然R的值增大了，但是G/B值的大小也会影响混合的颜色，导致变成了灰色。所以RGB三个值，牵一发而动全身，如果想在RGB空间找到红色范围是非常困难的，这就需要将色彩从RGB空间转换到HSV空间（如下图，这里不做详述），在HSV空间红色的范围可以很容易的表示出来。RGB空间就和时域一样，都有着自身的限制。所以最容易理解的表现形式并不一定是最方便计算的。我们往往需要进行一种变换，将在原来空间中难以处理的问题变换到方便计算的空间中去。频域频域就是描述频率所用到的空间或者说坐标系。频率虽然比较抽象，但是在我们的生活中是无处不在的，只是我们很少直接提到这个专业名词。对于波来说，频率是每秒波形重复的数量。声音是一种波；光具有波粒二象性，也具有电磁波的性质；更普遍的说，频率是物质每秒钟完成周期性变化的次数。比如家里用的交流电是50Hz，意思就是电压每秒完成50次振荡周期，如下图。交流电波形而前面提到的低音效果是什么样的效果呢？就好比家庭影院中的低音炮，它是如何实现重低音的呢？简单来说，可以将它简化成一个低通滤波器，下图是低通滤波器的频率响应曲线。低通滤波器频响曲线横轴是频率(Hz)，纵轴是声音大小(dB)。（请忽略图中的频率刻度，没有对应人声的频率范围）所谓的低音效果，其实就是对人声中的低音部分保留或增强，对应上图中左侧的横线部分；而对于人声中的高音部分进行衰减，对应上图中右侧的斜坡部分。通过这个低通滤波器，我们就能将低音过滤，将高音衰减。为了实现更好的视听效果，实际中，功放或播放器的实现会比这个复杂得多，上图中进行了极简化。可见，低音效果是在频率范围内考虑问题，而波形图是在时域内的图像，所以如果想在时域内解决低音效果的问题，就如同鸡同鸭讲。所以我们要就要找到一个沟通时域和频域的桥梁，也就是一个翻译，让时域和频域能够无障碍的沟通。但是，时域和频域表达的又只能是同一种信息，只是表现形式不同。就好比人们想了解古埃及文化，但完全不了解古埃及象形文的含义，所以也就无法根据记载的文字了解当时的文化。直到商博良破译了罗塞塔石碑上的古埃及象形文，才打开了古埃及文化的大门。所谓破译，其实就是找到古埃及象形文是如何表意的，然后翻译成现有文字系统，比如希腊文。它们本质上表达的是同一种信息，只是表现形式不同。时域转频域极坐标与直角坐标系类比前面类比了RGB空间，解释了为什么要进行时域到频域的转换。可能还不够形象，这里再用直角坐标系和极坐标系做一个类比。我们来看一下阿基米德螺线(如下图)，当一点P沿动射线OP以等速率向外运动的同时，这射线又以等角速度绕点O旋转，点P的轨迹称为“阿基米德螺线”。它的极坐标方程为: r = a + b\theta 。这种螺线的每条臂的间距永远相等于 2\pi b 。这种曲线在极坐标系中很容易的表示出来，而且形式非常简单优雅。但是在直角坐标系下要以X-Y的形式表示出来确是非常困难的，只能用参数化方程来表示。也就是说，有些问题，当我们换一个空间或者说域去考虑的时候，可能会豁然开朗。傅里叶级数为了形象的理解为什么要进行时域到频域的转换，前面已经举了很多的例子，下面正式开始进入时域和频域的变换。我们先来看一下标准正弦函数，如下图。在时域它的函数方程是 y = \sin(x) ，而它的频率是 f =\frac{1}{T}= \frac{1}{2 \pi} 。所以，上面这个函数在频域中的图像如下频域函数横轴是频率f，纵轴是幅值A。上面两张图分别从时域和频域展示了正弦函数，但表达的都是同样的信息。更一般的有 y = A\sin(2 \pi f x + \phi) , 其中 f
是正弦函数的频率，\phi 是初始相位， A 是幅度。在广义的频率中， f
可正可负，上图中旋转臂顺时针旋转， f 为负值。如果旋转臂转的越快，则频率越高；零时刻旋转臂和水平方向的夹角，就是初始相位。由于正弦函数是单一频率，在频域中只需要一根竖线就能表现出来。我们期望的也是将时域的信号转换成一个个单一频率的正弦函数的组合，这样我们就能够在频域中用一根根竖线表示出来，也就完成了从时域到频域的转换。而上面提到的正弦函数表达式可以转换成如下形式，{\begin{aligned} y &= A\sin(2\pi f x + \theta) = A\sin(\theta)\cos(2 \pi fx) + A\cos(\theta)\sin(2 \pi fx) \\&= a_n\cos(2 \pi fx) + b_n\sin(2 \pi fx) \end {aligned} } 所以，如果可以将任意波形都转化成若干个正弦函数和余弦函数的线性组合，我们是不是就完成了时域到频域的转换？别急... 实际的波形可能会有一个”直流分量“，如下图。这个方形波并没有沿X轴往复运动，而是沿 Y=2 这条直线往复运动。对于这类的波形，单纯的用正余弦函数组合是无法表示出来的，因为正余弦都是沿X轴往复运动，所以必须叠加一个”直流分量“。所以最后，如果任意波形都可以转化成常数、若干个正余弦函数的线性组合，我们就可以完成时域到频域的转换。用数学公式表达如下面所示：{\begin{aligned}s_{N}(x)& \stackrel{?}{=}
{\frac
{a_{0}}{2}}+\sum _{{n=1}}^{N}\left(\overbrace {a_{n}}^{{A_{n}\sin(\phi _{n})}}\cos({
{2\pi f nx}})+\overbrace {b_{n}}^{{A_{n}\cos(\phi _{n})}}\sin({
{2\pi f nx}})\right)\\&\end{aligned}} 上式中的 \frac
{a_{0}}{2} 就对应了直流分量，我们可以把它想象成一个常数而已。于是问题就转化成，对于任意波形，我们能不能找到一组系数 a_n 和 b_n ，使上述等式成立？（为什么上式采用了离散的频率，而且都是 2 \pi f 的整数倍呢？后面会介绍）。到这里，法国数学家傅里叶就必须登场了。他在1807年发表的论文中帮我们完成了这个工作，他提出了一个当时非常具有争议性的论断：任何连续周期信号都可以由一组适当的正弦曲线组合而成。其实，对于连续周期信号，比如上图中的周期方波，严格意义上说它的频域变换叫做傅里叶级数，因为经过频域变换后，它的频谱是离散的。而当我们现在说起傅里叶变换，默认指的是连续非周期信号的变换，如下图所示。因为非周期信号可以想象成信号的周期趋近于无穷大，所以傅里叶变换其实是对傅里叶级数的扩展。非周期信号正交性我们接下来介绍的都是基于连续周期函数的频域变换，也就是傅里叶级数。重新复制一下前面要证明的等式。{\begin{aligned}s_{N}(x)& \stackrel{?}{=}
{\frac
{a_{0}}{2}}+\sum _{{n=1}}^{N}\left(\overbrace {a_{n}}^{{A_{n}\sin(\phi _{n})}}\cos({
{2\pi f nx}})+\overbrace {b_{n}}^{{A_{n}\cos(\phi _{n})}}\sin({
{2\pi f nx}})\right)\\&\end{aligned}} 在我们尝试求解系数 a_n 和 b_n 之前，我们先解释一下前面留下的问题。上式为什么采用了离散的频率，而且都是 2 \pi f 的整数倍呢？这样的一组正余弦函数除了可以表示单一频率之外，方便的在频域表示，而且组成一组正交基，具有两两正交的优质特性，可以方便的计算系数。我们知道，正交是线性代数里的概念，是垂直这一直观概念的推广。比如，在欧几里得空间中，正交就是两个向量的内积为零。如下式。{\displaystyle \langle (x_{1},\ldots ,x_{n}),(y_{1},\ldots ,y_{n})\rangle =\sum _{i=1}^{n}x_{i}y_{i}=x_{1}y_{1}+\cdots +x_{n}y_{n} = 0} 下图中，X、Y、Z三个轴就是两两正交的。因为X、Y、Z三个轴对应的向量分别是(1, 0, 0), (0, 1, 0), (0, 0, 1)，根据上面的内积公式可得他们两两内积为零。那对于两个连续函数来说，应该如何表示正交呢？函数在某个区间内部有无穷多个点，无法直接套用内积公式。但我们可以借鉴积分的思想，将函数在一段连续区间分割成一份一份，这样每一份的取值合起来就可以组成一个向量，于是可以用向量的内积来表示两个函数是否正交。如下图当分割的区间无限小时，向量变成无限维，于是向量的内积就可以用积分来替代了。所以两个函数的正交其实可以用积分来表示。对于 \sin4x 和 \sin2x ，在一个周期内，处于X轴上方的面积和X轴下方的面积相等，所以这两个函数的积分为0，也就是互相正交的。还有另外一种情况积分的两个函数相同，都是 \sin4x ，这时积分结果都在X轴上方，积分大于0，也就是不相互正交。更一般的情况，有\begin{aligned}
& \left\{
\begin{aligned}
&\int_{-\pi}^{\pi} \cos(mx)\, \cos(nx)\, dx = \pi, \quad m = n,
m, n \ge 1 \\ &\int_{-\pi}^{\pi} \cos(mx)\, \cos(nx)\, dx = 0, \quad m \neq n, m, n \ge 1 \end{aligned} \right. \\\\& \left\{
\begin{aligned}
&\int_{-\pi}^{\pi} \sin(mx)\, \sin(nx)\, dx = \pi, \quad m = n,
m, n \ge 1 \\ &\int_{-\pi}^{\pi} \sin(mx)\, \sin(nx)\, dx = 0, \quad m \neq n, m, n \ge 1 \end{aligned} \right. \\\\ & \int_{-\pi}^{\pi} \cos(mx)\, \sin(nx)\, dx = 0, \quad m, n \ge 1 \\& \int_{-\pi}^{\pi} \cos(nx)\, dx = 0, \quad n \ge 1 \\& \int_{-\pi}^{\pi}
\sin(nx)\, dx = 0, \quad n \ge 1 \end{aligned} 所以， 1, \cos(x), \sin(x), \cos(2x),\sin(2x), ... , \cos(nx), \sin(nx)
不仅能够表示单一频率，而且还构成了一组正交基。于是我们重新看上面要证明的等式，如下所示{\begin{aligned}s_{N}(x)& \stackrel{?}{=}
{\frac
{a_{0}}{2}}+\sum _{{n=1}}^{N}\left(\overbrace {a_{n}}^{{A_{n}\sin(\phi _{n})}}\cos({
{2\pi f nx}})+\overbrace {b_{n}}^{{A_{n}\cos(\phi _{n})}}\sin({
{2\pi f nx}})\right)\\&\end{aligned}} 上面这个等式和前面提到的正交基不同的仅仅是频率 2 \pi f ，当积分区间修改成对应的周期范围内，则所有的特性完全一样。也就是说如果上面这个等式能成立，将是一个非常完美的等式。我们不仅完美的将它转换成了频域，用单一频率的正余弦函数组合表示出来，而且他们还是一组正交基，从线性代数的角度来看就是线性无关的，彼此互不影响，就和欧几里得空间中的X-Y-Z轴一样，无论我如何修改X的值，对Y和Z都是没有任何影响的。这种正交的特性可以让我们非常方便的求出对应的系数。比如我想求出 a_{10} ，它的系数就是 \cos(2\pi f \times 10 \times x) ，前面已经提到， \cos(2\pi f \times 10 \times x) 是正交基中的其中一个基底，它只与自身积分不等于零，而与正交基中任意其它基底积分都为0。有了这个非常好的特性以后，我们要求 a_{n} ，只需在两边同时乘以 \cos(2\pi f \times 10 \times x) ，然后做积分，其它所有的频率部分因为正交性，都变为零，等号右侧只保留了 a_n 的部分，我们就可以求出 a_{10} 。更一般地，对于 a_n 和 b_n ，有如下等式成立{a_{n}={{2}{f}}\int _{x_{0}}^{x_{0}+\tfrac{1}{f}}s(x)\cdot \cos({ {2\pi f nx}})\ dx} \\ b_{n}={{2}{f}}\int _{{x_{0}}}^{{x_{0}+\tfrac{1}{f}}}s(x)\cdot \sin({{2\pi f nx}})\ dx \\ 到这里，我们利用正交性求出了傅里叶级数中的 a_n 和 b_{n} ，求解过程并没有那么严谨，只是能够直观的理解如何进行的时域到频域的转换，以及如何利用正交的特性去求解系数。过冲现象虽然求解出来了，但是真的如傅里叶所说，我们可以用正弦波去表示任意的连续周期函数吗？以方波为例，我们看一下下面的动图。方波如上图，随着频率越来越丰富，合成的波形也越来越接近方波了，当n趋近于无穷大，也就是频谱范围无限大的时候，就可以无限逼近方波了。但是，我们注意到，即便在n = 29的时候，合成的方波还是棱角处还是有一定的过冲。而且我们发现，在n增大到29的过程中，这个过冲并没有明显的减小。那么在n趋近于无穷的时候真的能够避免这种过冲吗？我们知道，正弦函数是一个处处连续且可导的函数，也就是说正弦函数是一个比较圆润的函数；而方波却是有棱角的，在棱角处是不连续的。一个圆润的函数最后可以合成一个带有棱角的函数吗？即使将频谱范围扩展到无穷大，就真的能够逼近出棱角吗？事实上，当初拉格朗日也是拿这一点反驳傅里叶的。在不连续点的过冲即使频域扩展到无穷，过冲也不会降为零的。其实傅里叶所说的逼近其实是能量的无限逼近，也就是经过傅里叶变换后的波形能量和原始波形能量可以无限逼近。复频域傅里叶级数到此，基于三角函数形式的傅里叶变换已经介绍完毕。但傅里叶变换还有一种复频域的表示方式，通过复频域表示更加简单直观，但这就需要用到大名鼎鼎的欧拉公式。对任意实数x，都存在
e^{ix}=\cos x+i\sin x ，
i 是虚数单位通过复频域表示时，会出现虚数单位，而这个是我们在三角级数表现方式中不曾出现的，而最后复频域表示方式要能够化成和三角级数的相等的表达形式，所以必须想办法消掉虚数单位。所以我们就想到共轭复数
e^{-ix}=\cos x-i\sin x 。有了共轭复数，我们可以通过两个互为共轭的复数加法将虚数消掉。于是有下面的式子，我们将频域的 1 ~ N求和 和常数 转化成复频域的 -N ~ +N求和，这样通过构造 -N，就会出现两个互为共轭的复数。 {\displaystyle {\begin{aligned}s_{N}(x)&={\frac {a_{0}}{2}}+\sum _{n=1}^{N}\left(\overbrace {a_{n}} ^{A_{n}\sin(\phi _{n})}\cos({{2\pi f nx}})+\overbrace {b_{n}} ^{A_{n}\cos(\phi _{n})}\sin({{2\pi f nx}})\right)\\&=\sum _{n=-N}^{N}c_{n}\cdot e^{i{{2\pi f nx}}} \\&
=\sum _{n=-N}^{N}c_{n}\cdot (\cos({{2\pi f nx}}) + i\sin({{2\pi f nx}})) \end{aligned}}} 然后我们用待定系数的方式求解 c_n ，另 c_n = \left\{ \begin{aligned}
p + i q \quad (n>0) \\ p-iq \quad(n < 0) \end{aligned} \right. ， 将互为共轭的两个系数提取出来，对比三角级数和复频域的表示方式，列出等式。\begin{aligned} &(p - iq)(\cos({{2\pi f nx}}) - i\sin({{2\pi f nx}})) + (p + iq)(\cos({{2\pi f nx}}) + i\sin({{2\pi f nx}})) \\=& a_n\cos({{2\pi f nx}})+b_n \sin({{2\pi f nx}})) \end{aligned} 解得 p = \tfrac{a_n}{2}, \quad q = -\tfrac{b_n}{2} ，于是有如下表达式c_{n}\ {\stackrel
{{\mathrm
{def}}}{=}}\ {\begin{cases}{\frac
{1}{2}}(a_{n}-ib_{n})&{\text{for }}n>0\\{\frac
{1}{2}}a_{0}&{\text{for }}n=0\\c_{{|n|}}^{*}&{\text{for }}n<0.\end{cases}} 所以我们有了复频域的傅里叶级数表示方式。如果想要求出 c_n ，同三角级数一样，在复频域上 e^{i{{2\pi f nx}}} 同样具有正交性，所以我们想要求出 c_n ，只需要在等式两边同时乘以 e^{-i{{2\pi f nx}}} ，然后再进行积分，就可以过滤掉其它复频率分量，而只保留 c_n ，于是我们有c_{n}={f}\int _{{x_{0}}}^{{x_{0}+\tfrac{1}{f}}}s(x)\cdot e^{{-i{
{2\pi f nx}}}}\ dx结尾傅里叶就是在它的《热的解析理论》中提出了傅里叶变换的一系列思想，虽然他如此伟大，但是他最后的结局却是“no zuo no die”排行榜第一。傅里叶对热极为痴迷，同时认为热是世界上的最棒的东西，甚至可以包治百病！他为了证明这个理论，一次他在身体不舒服的时候，在大热天，他把门窗四闭，烤着火炉，“治疗”着自己。变态的室温大大加重了他的病情，最终活活自己热死了。。。1830年5月16日，傅里叶卒于法国巴黎。欣赏到此，我们已经介绍完了傅里叶级数的三角函数和复频域形式。最后，我们欣赏一下傅里叶变换都能够模拟什么样的波形。如果你觉得傅里叶级数只能画上面的“心电图”，那你太小看它了，我甚至可以用它来画恐龙，下面这张图是由50个频率的傅里叶级数组成，在我本地跑了将近4分钟（实在太卡），20倍速播放。恐龙画一只可爱的小猫也不是不可以~小猫咪傅里叶发现了这么伟大的公式，连接了频域和时域，其实给我们带来的更多的是思维上的转变。当我们在时间轴上思考问题比较困难时，是不是可以换一个维度重新审视思考或许会豁然开朗。}

序言每增加一个数学公式都会使读者减半，原话出自霍金的时间简史。所以我知道这篇可能不会有多少人看，但是仍然执着的想分享出来。（长文多图，流量慎点）佛语中有箴言：坐亦禅,行亦禅,一花一世界,一叶一如来世间万物是复杂的，但是又是纯粹的简单的，从宏观的花花世界到微观的原子电子，万物都在按照它的规律运行，而我们的先辈前人，一直都在用自己的方式与经验，总结着万物运行的规律。不要误会，本文并不是哲学论题的讨论，一个偶然的机会拜读了知乎上的一篇关于频域水印的问答。（阿里巴巴公司根据截图查到泄露信息的具体员工的技术是什么？ - 知乎）当中有对频域数字水印的实现与讨论，身边有不少的朋友对此颇感兴趣，于是我就想以后机会写一个“从零开始的频域变换到水印的完整解答”，当然，你也可以直接跳转到本文的最后一章节“鲁棒盲水印”来查看频域水印是如何在版权保护中起作用并对抗各类版权偷盗者的攻击的。本人并非全才，本文中多有疏漏错误还望各位读者多多指正，文章将从最基础的三角函数开始，一步一步地推到并演化到傅里叶变换直到频域签名。我相信学习如流水，希望读者能在本文的循序渐进的推导过程中满足对相关技术的求知欲，并对文中的不足与错误不吝赐教。最后再次引用箴言的后半句作为导言的结束语：春来花自青,秋至叶飘零,无穷般若心自在,语默动静体自然。为了避免一开始就引入那些过于严肃的话题，笔者希望在开篇的部分，做个简短的说明来告诉大家，这篇文章做的是什么的。首先，笔者画了一张图名叫：《虎虎生威》，好吧，就是下面这张因为笔者太叼了，所以画的图应该署名一下，你可以看到右上角的DBinary，没错，那就是笔者的大名，但是，盗图狗很快就把我这张大作给偷走了，最后他居然把自己的名字写了上去好气啊，明明是笔者画的图，怎么变成张三了，现在死无对证了，到底是谁画的，笔者暗暗不爽，于是做了一台时光机，回到笔者画图后还没发布之前，不行，这回不能光靠签名，要加点靠谱的水印，于是，笔者开发了一款软件，对这个图片进行了隐签名：（这是签名后的图片，好吧，因为颜色过于单调还是能看出明显干扰的，实际对照片签名几乎看不出来）盗图狗果然还是出手了，它篡改了我的图片我一看跳了起来，***张三你怎么盗我图，显然，张三不服，凭什么说我盗你图，证据呢？？？我不慌不忙打开频域程序，加载了张三盗窃后的图那么张三同志，麻烦你告诉我图像频域里为什么写的是我的名字？张三哑口无言.....这就是本文将要讨论的技术细节，如果阅读到这里你还不明白本文的主旨与目标是什么，你可以直接跳转到本文的最后一章节“鲁棒盲水印”来查看频域水印是如何在版权保护中起作用并对抗各类版权偷盗者的攻击的。从三角形开始有人说，上帝使用三角形创造了这个世界，这点我完全同意，三角形拥有如此之多的特性，足够让每一个探究者为之所着迷，它仿佛维系着几何与世界的基石，诞生出数学中众多的定律并在今天成就了我们的学术与技术大厦。当然，本文并不是抒情散文，我并不打算也没有这个能力去探究那些更深层次的数学理论关系，但理解三角形并不是制作变形金刚，你可以在一张纸上画三个点然后用直线把他们连起来，那就是一个三角形了如图（a.1）三角形有非常多的特性，首先，确定它是一个稳定的结构。另外确定一个平面也仅仅只需要一个三角形就足够了，三角形的所有内角角度之和是180°（a.2）。但最为有意思的是被称之为直角三角形的东西，在直角三角形中有一个角的角度是90°（a.3），例如图（a.3）就是一个直角三角形。实际上在这个直角三角形中，三角形的边长比例，将随着角度θ的变化而变化，一个角度θ决定了abc三边长度的比例关系，于是在这里，我们引入了一个叫三角函数的东西，其中，，，当θ小于90°的时候，我们可以在直角三角形中非常直观地看到三角函数的变化关系，但在直角三角形中，θ的取值范围，也被限制到了0到90°，为了能让θ表示更大的范围，我们就需要引入新的表达方式了。我们首先先建立一个直角坐标系（a.4），然后以原点为圆心，画一个圆。同时，我们在圆上任意取一点，并将该点的坐标设为（x，y）通过勾股定理我们知道，圆上的点到原点的距离，r=那么，对三角函数我们重新定义为，，，这样一来，我们就可以表示θ为任意实数时三角函数对应的值了通过这个图我们同时可以知道，当时，实际上相当于在一个圆上绕了若干个个圈，你可以想象看着你家里的时钟，秒针每过60秒就转了一个圈，你现在看到秒针的位置，在60秒后，120秒后，180秒后…它仍然会指向同样的位置，这个特性在三角函数中同样有效，我们管它叫做三角函数的周期性。现在，让我们引入一个新的符号π，π在角度上代表180°，除了角度外，我们还引入弧度，它的定义是弧长等于半径的弧，其所对的圆心角为1弧度。实际上半圆弧长和半径的比例恰好是一个定值，因此，π除了在角度上表示180°，在弧度上则是一个定值，这个值是一个无限不循环小数，我们常常将它约等于为3.14。所以，同时因为周期性，其中n是一个整数。那么，在第一章节，三角函数的几何意义，就显而易见了。让三角函数动起来假如我们把时间引入进来，那么，三角函数就开始变得生动了，最直观的比喻就是现在挂在大厅墙上的时钟，秒针每分钟都会转动360°，现在，让我们假设秒针指向“12”，也就是垂直于水平面时，它的角度是0°，那么经过15秒后，秒针指向“3”，也就是转动了90°，30秒后经过了半分钟，指向“6”也就是转动了180°，那么，秒针每秒转动的角度是,我们用t来表示时间，那么时间与秒针转动角度我们可以用θ=6t来表示，在物理上，我们常常使用ω来表示角速度，那么，时间内转过的角度就是θ=ωt现在，让我们画一个二维坐标系，并设横坐标为时间，角速度ω=2π，那么，的坐标如下图（b.1）所示可以看到，在一秒0-1秒，已经是一个完整的周期了，因为它在1秒内拥有一个完整的周期，我们就称他的频率为1hz，显然的当ω=4π时1秒内有两个完整的周期，因此，它的频率就为2hz（图b.2）。可以看到，频率实际上和角速度是相关联的，我们使用字母f来表示频率，函数公式如下角速度和频率是一个正比关系，角速度越大，频率也就越大对正余弦函数而言，我们也常常写作三角函数的一些常用公式我们可以非常直观的从几何图像中推导出三角函数的一些性质，画一个坐标系，同时以原点为圆心绘制一个半径为1的圆（图c.1）那么，可知由周期性可知同时，三角函数间是可以互相转换的，我们很容易得出这样的公式：因为正余弦函数存在这种互相转换的关系，因此之后我们统称它们为正弦函数我们还可以进一步证明二角的更多公式（证明过程略），这些公式我们都可以直接拿来使用二角和差：和差化积通过联立二角和差公式，我们还可以得到积化和差公式和二角和差推导出的二倍角公式三角函数的正交性在讨论正交性之前，让我们用最简单的方式了解下微积分，我们取一个最简单的一元函数做比方现在，让我们画一个坐标轴，那么，这个函数的图像是这样的(d.1)：现在，做一条经过（2,0）并垂直于x轴的直线，交于点A，（2,0）为B，原点为C如图（d.2）那么，我们很容易求出三角形ABC的面积实际上这个求面积的过程，就是微积分于该函数上的几何意义，这个求面积的过程，我们用微积分来表示，就是同样类比的，这个积分方程实际上是求四边形ABCD的面积（d.3）在二维平面上我们很容易将微积分理解为函数在一定范围内和x轴围成的面积，在但实际上面积和微积分稍有不同，因为面积肯定是一个非负数，但积分却可以是负数，例如他在坐标系中是三角形HCI的面积，但是它的x轴下方，我们可以理解为三角形的高度是一个“负数”，因此，这个区域的面积也是一个“负的面积”，所以，这段的积分为-2（d.3）那么，假如我们对这个函数的-2到2积分，那么就会变成一个正的面积加上一个负的面积，结果它们相互抵消了，结果变成了0，如下推导实际上连续的中心对称函数图形h(x)，在的积分都是0显然的，正弦函数图像满足这个性质（d.4），不仅如此，我们可以非常直观的看出来，正余弦函数在函数的一个周期内的积分都是0（d.4 d.5）现在，我们假设有两个自然数m，n且，并假设两个函数然后我们将这两个正余弦函数任意进行组合，并对他们在-π到π进行积分通过上述的常用公式，我们可以推导出下面的三个式子可以看到，m不等于n时，积分结果都是0从正交性得到的启发还记得之前我们使用的将时间作为变量的三角函数么现在，因为三角函数的特性，我们管三角函数在时域上的表示叫波，当然也就是正弦波了，我们用更加通用的公式来表示一个正弦波\或者，余弦波也同样是这个道理，毕竟它与正弦波的不同仅仅只是“移动了一下”我们很容易证明这点。现在我们可以想象一个正弦波在二维坐标上的样子，或者是我们在科幻电影或科研实验室中，看到仪器仪表上那一段段的波形。或者更加直观的，我们将一块石投进水池里，荡起的波浪也像极了正弦波。当然，在自然界一般不会出现标准的正弦波，各种波的叠加，你可以想象在一个房间里一群朋友聊天，或者是在KTV中歌声和谈话声混杂的场景，它就是声波的各种叠加尽管聊天中大家都在说话，但是我们仍然不会把朋友们说的话混淆起来，因为不同的人发出的声音频率也不一样，这也就是我们常说的未见其人先闻其声，我们天生具备有分别不同频率声波的能力，所以尽管环境嘈杂，我们仍然能够取得我们想要的信息。但是现在我们如何使用数学的手段，查看波形函数f(x)中是否有我们想要频率的波形呢。显而易见的，答案当然就是本章节的标题，应用波的相关性，例如我们使用一个正弦波去乘以一个由多个正余弦波叠加而成的函数f(t)，然后对它们进行积分从相关性可以知道，当角速度（上一章节的m，n）不同时，也就是频率不同时，积分结果是0，只有频率相同时，积分结果才不为0，因此，如果f(t)中包含1HZ的正弦波，那么积分的结果就不为0，否者，结果就是0那么，一个检波手段也就诞生了。欧拉公式复变函数正弦函数和复数表面上并没有什么关联，但正所谓千里姻缘一线牵，在我们考虑着那根看不见的红线另一端系的是谁的时候，欧拉早在几个世纪前就将正弦函数与复数牵了一条红线，从而撮合了数学上这一段流传千古的因缘，不过如果再在红线上扯下去，这篇文章就要变爱情小说了，但是我们仍然需要回到这个渣作的现实三次元，现在让我们来介绍一下大名鼎鼎的一个复变函数，它的知名程度基本上在是数学上的安徒生童话，人人皆知它的公式如下其中，i表示复数的虚部，它是，也就是说，当然我们并不能深究他在现实生活中的意义，但它在数学多个方面，起着举足轻重的意义。那么，复数如何理解呢我们来看看复数的标准形式其中，a，b为任意实数，a在复数中表示实部，b在复数中表示虚部假设我们现在画一条二维坐标系，x轴为实部，y表示虚部，那么，1+2i实际上表示的是坐标上的（1,2）那么在坐标系上，实际上是一个半径为1的圆（d.7）复数在信号的表示现在让我们来考虑下面这个复数它的实部为，如果我们希望用的形式来表示它，那么，它就变成了画在坐标系中，如图d.8按照极坐标的角度而言，即一个长度为2的变绕X轴逆时针旋转现在让我们更进一步放大我们的脑洞，以原点为圆心，2为半径画一个圆（图d.9）那么,假如我们将看作是与原点距离为1，绕x轴逆时针旋转的点，将与之相乘就变成了，初始位置为(以角速度以原点为圆心逆时针旋转的点那么，假如我们将复数用在信号中，就从复数坐标系上的点拓展到了在线信号的幅度与初相了。还记得之前我们写的通用的正弦/余弦函数么这就是上面的通用公式，因为这个公式太重要了（或者说避免你翻回去找这个公式），我将它再写了一遍可以看到，对于一个正弦函数，我们仅需要知道它的幅度，角频率，初相，就可以确定这个正弦函数是什么样子的了但正弦函数比善变的女人还有能耐，通过一些数学变形，你还会发现它有时比哄妹纸有意思多了，通过三角函数的公式，我们可以得到那么，这个函数就变成了正余弦函数的组合，同时我们也得出了那么就有傅里叶的故事这里我们抛开那些繁琐的数学公式，然后把时间拉回到17世纪，当时有一位法国的男爵名叫巴普蒂斯·约瑟夫·傅里叶（Baron Jean Baptiste Joseph Fourier）当然，并不是因为它当了一官半职我们才介绍它，但他的成就，却从17世纪一直影响到我们今天，毫不夸张的说他奠定了信号系统的基础(或者说给出了指导方向？)。具体时间还是要回到18世纪，伯努利（D.Bernoulli）就曾经提出：一个弦的实际运动都可以用标准张模的线性组合来表示，但是当时另一个数学大神拉格朗日是不兹磁这个说法的，拉格朗日认为，不可能用三角级数来表示一个具有间断点的函数，就在这个环境下，傅里叶仍然坚信：“任何”周期信号都能够成谐波关系的正弦级数来表示，他还专门写了一篇论文，但迫于拉格朗日当时在学术界的威望，傅里叶理论的论文一直没能被发表，直到傅里叶的晚年，他才得到他应有的承认。当然，以我们现在的角度来说，傅里叶当时的理论是有缺陷的，在后人不懈努力下，傅里叶变换才趋于完善，实际上傅里叶变换并不是因为傅里叶对这个理论在数学上做了多大的贡献，但傅里叶当时对问题的前瞻性与指导性，却奠定了这个数字信号举足轻重的公式基础。检波与傅里叶变换现在我们讨论的就是傅里叶大神“任何连续周期信号都可以由一组适当的正弦曲线组合而成”这一命题，在信号系统或者是工学，傅里叶变换绝对是入门的加减乘除，实际上傅里叶变换并不是一个公式，而是多个分别应对不同类型信号的多个公式，但万变不离其宗，不管变体如何变换，其核心的思想是不会改变的。那么傅里叶变换的意义是什么呢，回想一下我们之前说的正弦函数的相关性，通过这个相关性，我们可以检测某个波里是否包含某个频率的正弦波，进一步的，假设我们用无限多个不同频率的正弦波与它正交组合，我们就能知道，这个波是由哪些频率的正弦波组合而成的了，我们管波在时间轴上的表达，叫做时域，图b.1其实就是的时域图，通过傅里叶变换，我们能够将是由哪些频率组成的图像表示出来（当然只有一个1HZ频率），也就是横坐标由时间变成了频率，也就是我们说的频域了。简单来说，在时间信号处理中傅里叶变换就是将时域信号转换为频域信号的一个变换过程这里，我们先来看看连续傅里叶变换，因为傅里叶英文的开头是F，所以我们使用大写的F来表示傅里叶变换，f(t)用来表示某连续非周期性的时域信号函数看看我们的欧拉公式然后把欧拉公式代入傅里叶变换角速度乘以时间，不就是了么，我们将这个变换函数写成这种形式再将带入公式得到结果变得显而易见了，简单来看，傅里叶变换与检波的手段多多少少的相似性---对sin与cos分别检波（为什么要使用同一频率sin，cos进行检波呢，因为对于一个频率的正弦波我们不仅仅要知道其存不存在，更需要知道其波幅及相位，用两个不同相位的正弦波进行检波，就可以取得其相位，在之后章节会进一步讨论）傅里叶级数的三角函数表达但仅从上面的检波手法来推断傅里叶变换，是不严谨的。现在设一个函数由一个直流分量（简单来说就是一个常数）和多个正弦函数组成，那么它可以写成这种形式int sum;for(n=0;n<N-1;n++){Sum+=n;}由上式可知，表示这个三角函数的角速度或者叫角频率，当n=1时，我们管叫基频，管叫傅里叶级数（余弦信号形式）利用三角函数的变换公式，上式可变形为设那么，上式变为现在，让我们正式的引入正交性的性质，还记得检波手段么，这里，我们假设对f(x)用sin(nwt)进行检波，那么就有假设f(x)中含有角频率的正弦波系数为，那么根据三角函数的正交性，上式就有进一步计算，可得周期连续时间傅里叶级数现在，让我们来想象一个函数f(x)它是一个周期函数，那么根据傅里叶的理论，它能够表示成若干个（无穷个）一组“适当”的正弦曲线组合而成，在前面几个章节，我们通过欧拉公式，得到显然的，这个复指数信号的频率是，现在我们假设有另一组“适当”的正弦函数，它们的频率刚好是的整数倍并且幅度也不同，那么，这组信号我们可以使用来表示（k为自然数）。那么从上面的根据欧拉公式，我们也很容易得出下面的推导现在将这两个公式带入得到化简后得因此，上式最终变为了设上式就写成了,n=k并对两边同时在一个周期内积分，那么我们就得到公式\进一步变形为于是，得到也就是实际上an就是我们所说的傅里叶级数，或者说是频域系数。通过这个系数的值，我们可以知道这个频率的波对原始波的能量贡献值，在之前的《信号的复数表示》章节中，我们可以了解到这个系数确定了该频率波的幅度，初相，从而完成信号时域到频域的分解，并且我们还知道了也就是说，通过ak的模我们知道对应角频率波的波幅（等于该频率幅值的一半）通过可以得出对应角频率正弦波的相角。周期离散时间傅里叶级数在上一个章节中，讨论了周期连续时间信号的傅里叶级数求解方式，那么，连续信号可以求其级数，离散的是否也有这样一个公式呢但在介绍离散变换变换前，我们先来了解连续和离散是什么，。其实顾名思义。比如给你一段长度100米的绳子，当然，这个100米的绳子是连续的，如果你在50米的地方剪短这根绳子，那么它就是不连续的了，假设我们把这根绳子切成若干个段。直到每个段都变成一个“点”，我们可以直接用数字编号每一个点，那么他就变成离散的了。用图像来继续说明，如图e.1，这是一个sinx的函数图像，当然，它是周期无限长且连续的现在，我们每隔二分之π就取一个点，那么这些点在坐标系上就是周期离散的（e.2）说到这里，我们现在要介绍的变换，就是处理这些点的变换方式首先，因为是等间距对周期信号采样，所以采样出的点也是周期性的，假设采样的点用数组x[n]来表示，也就是说假设周期为N，x[a]与x[a+Nk]是完全相等的值，也就是说，整个离散样本中取任意周期N内他们的累加和是一样的，同样的，x[n]仍然可以用“恰当”的正弦函数组合而成（或者说可以由基波频率为的一系列波形组合而成）基于这点，我们就可以将x[n]写成这种形式（k=<N>表示取任意连续N点即一个周期内点，不管如何取，结果都是一样的）我们取可以知道，当K不等于N时的结果为0，仅当K等于N时的结果为N。同样的，我们再次将两边同时乘以得到然后再同时对两边进行N项上求和显然的，仅当k=r或n为N的整数倍时，右式才不为0，那么，上式就变为那么频谱系数求得！。非周期离散时间傅里叶变换如果严格来说，自然界大多是没有理想的周期信号的，那么，是否有办法处理非周期离散时间的信号么。我们来看看这样一个离散信号（图e.3），它只有一个脉冲取样这个脉冲信号仅在-3，-2，-1上是有值的，其余的值都是0那么我们是否可以把它当成一个周期无穷大的周期信号呢。我们先来看看上一节中周期离散时间傅里叶级数的分析公式假如把这个分析公式套用在上述的信号中，那么因为仅在三点有值（其它都是0）并且周期是无穷大那么我们就得到了公式当然，它和下面这个式子是等价的这样我们就将公式推广到了更加通用的公式类型现在定义函数那么我们现在再将ar代入原式中，得到从式中看出，随着N趋近于无穷大，趋近于无穷小，那么，上式就从累加变成了积分，且因为的周期为2π，且其仅在周期内有值，于是，上式也随之变为了非周期离散有限长度傅里叶变换最后，我们将上述的变换公式进行进一步的推广，就是非周期离散有限长度傅里叶变换了，实际上它与周期离散傅里叶变换已经非常的接近了如果我们将有限长的信号推广到无限长的信号，那么我们先假设信号的样本点数为N个，那么，信号的n取值范围就可以定义在[0,N-1]我们假设将这个有限长的区间补到无限长，除了[0,N-1]区间，我们仅仅需要在其他区间再补上N个同样的离散信号就行了，这并不影响其结果，那么我们就可以将有限长非周期离散信号变为周期离散信号了，这样我们就可以直接套用周期离散时间傅里叶变换的分析公式：为了方便计算，我们周期取[0,N-1]，那么，公式就变成了当然在实际应用中，我们常常设：那么就有了有限长非周期傅里叶变换的的分析公式：非周期离散时间傅里叶变换的应用在上面几个章节中，我们从周期连续时间的傅里叶级数逐步的推广到了有限长离散时间傅里叶变换。虽然内容不少，但是真正实际在日常用的多的，是有限长非周期离散时间傅里叶变换，为何？当然我们幸运的不是在一个老牛拉破车的年代计算只能靠人脑算盘，现在的信号处理除非一些数学推导应用，大多数实际生活应用都是靠计算机来完成的，而这也决定了我们的公式必须与计算机的硬件相关，我们的计算机目前存储空间是有限的，而连续的信号（不管是频域还是时域），就相当于由无穷多个点组成，很遗憾，现在仍然没有存储容量无穷大的内存，就算有，也没有能够处理无穷大数据的CPU，因此，我们无法直接处理连续的信号，观察上面几个变换，也就只有时域和频域都是有限长度且离散的有限长离散时域傅里叶变换能够满足我们的需求了。下面为了说明方便，我们将有限长度且离散的有限长离散时域傅里叶变换的分析公式统一简称为离散傅里叶变换（Discrete Fourier Transform）或者其英文缩写DFT，将它的逆变换（Inverse Discrete Fourier Transform）也就是综合公式简称为IDFT。现在让我们来看看离散傅里叶变换对其中，N表示信号的采样点数，n的范围是0到N-1，的第一个点，的第二个点，以此类推，k的取值范围是0到N-1，由f=可知，其分辨率等于为了更加方便浅显地了解其中的计算，我们继续观察公式中的仍然是我们熟悉的味道，现在我们使用欧拉公式替换它于是离散傅里叶变换的公式变成了当k=0时可以看到，它实际上是所有样本点的累加和，那么它意味着什么呢，我想象一个波形函数它表示基频成分，因为我们从公式中可以看到，频率都是其整数倍，至于多少倍，就是由k值决定的。现在，让我们用一个实际的范例来验证离散傅里叶变换假设正弦函数我们假设对这个函数进行采样，采样频率是4HZ，那么，实际上我们将在0.25s，0.5s，0.75s，1s处取得其样本点，那么，对应的值应该是现在对其进行离散傅里叶变换那么因为有四个样本点，所以，N的值是4，k的取值范围是0-N-1，也就是0到3了我们先来计算的值这点没错，因为没有直流分量，所以它理所当然是0继续是的值也确实是1HZ。现在是的值这是2HZ的值，显然，它为0最后x3值是的值是2，但是我们知道，并不包含3HZ的波，实际上根据香农采样定律，4HZ的采样只能表达2HZ的波，因此这个点实际上是不准确的。但是，出现2的结果并不是偶然，我们接着往下看。巧妙的对称性，共轭如果说对称是世界上共有的一种美的表达，那么，在几何平面上，无穷多的函数就拥有这种对称性，在对称美这一个方面，出色的数学家也许并不逊色于毕加索。当然深究的话就是后话了，在这里我们来简单看看几个拥有对称性的简单函数：这是一个非常简单的二元一次方程，可以从图f.2中看到，他是关于y轴对称的，这句话如果用数学的语言来将，非常简单且直观的，我们可以用下面的公式来表达这个“对称”的思想对于这类关于满足上述公式的函数，我们管它叫偶函数。现在再让我们看另外一个函数它的函数图像如图f.3可以看到，函数的图像是原点对称的，相等，我们用如下公式表示这种“原点对称”关系我们管它叫奇对称当然，对称未必一定要是y轴或者是原点，现在我们将的函数图像向右平移两个单位，变成这样，它就关于x=2这条竖线对称了（f.4）这回我们在数学上用这个方程式来表示它关于x=2对称通过变形，它可以写成最后，我们用更加通用的方程式来表示某个函数f(x)关于x=N/2这条垂线对称现在，让我们来看看共轭是什么：你看过几米的《向左走，向右走》么，共轭是数学中一个文艺而浪漫的代名词，一句相当文艺诗来总结这个关系，就是：向左走，向右走，纵使背道而驰，相隔万里，但只要心连接在一起，也终会在大洋的彼岸，迎来相会的交点。文艺的诗歌艺术生看到后也许就开始感叹，多么美的诗句，世间万物芸芸众生，千里有缘一线牵，感谢我能遇见你，但对于地理大神也许不削一顾：这不就是说地球是圆的么，数学系的牛人毅然站出来，别做梦了少年，你们只会越离越远，就算到了世界末日，你们也碰不到一块儿。结果，文艺青年赢的了女神的芳心。理工大神则斩获了“屌丝”的称号实际上，共轭我们可以理解成相关联，也就是所谓的“缘分”，就是这一对数存在某种联系的意思，这里我们不深究更深层次的含义，我们主要说说共轭复数，首先，复数的形式是：它的共轭是非常简单的一句概括：实部相等，虚部取反，这两个复数互为共轭。那么，和我们之前说的对称性相结合的话什么是共轭对称性呢我们这样给出定义：当一个函数f其实部为偶函数，虚部为奇函数时，此函数就为共轭对称函数，即f(x)的共轭等于f(-x)，举一个非常简单的例子：显然的，实部是个偶函数，虚部xi是个奇函数，因此它是一个共轭对称函数，现在，我们来看看离散的范例，假设有以下复数序列我们分别提取实部与虚部那么分别是实部：1,2,3,3,2,1虚部：1,2,3，-3，-2，-1我们很清晰的看到了序列的对称性，实部偶对称，虚部奇对称。在这里，因为他是离散的序列，所以我们管它叫共轭对称序列。离散傅里叶变换的共轭对称性对称之美存在世间万物中的每一个角落，作为信号与系统中最美的变换函数之一，她也存在这种对称性。我们回到离散傅里叶变换的公式：我们假设：依据欧拉公式，可以得出那么，上式可以写成我们设k=N-k，并将它带入当中，那么就有和进行对比，发现，实部相同虚部相反，假如输入的信号为实信号时，刚好呈共轭对称性，将它代入离散傅里叶变换方程后其中*的意思就是共轭。也就是说离散傅里叶变换具有这种共轭对称性（输入为实信号时）。但共轭对称的范围是1~N-1，因此当k=0时，它的共轭对称并不存在（序列范围是0~N-1），所以我们需要额外讨论k=0时的情况：也就是之前所说的直流分量了。现在回到之前我们所求函数的序列离散傅里叶变换的结果根据共轭对称性得快速傅里叶变换在上文的范例中，我们仅仅是计算了序列长度只有4的离散傅里叶变换，也可以看到非常麻烦，倘若序列再长点的话，工作量就会呈指数级增长，我们使用矩阵运算来表达离散傅里叶变换的运算过程可以看到需要16次乘法运算，按照时间复杂度来算的话，它的复杂度是O()，当然，这还没算上sin cos与复数加减法带来的性能开销。所以，如果你想将一段2000长度的离散信号进行DFT运算，意味着你至少需要做400万次的运算，巨大的性能花销，导致了DFT在以前并不被看好，毕竟除非一些非常重要的信号，谁会花大量的人力物力去做这几百万次的运算呢，即使是在今天，几百万次的计算开销在有了计算机的帮助之后，也并不算一个小数目，倘若对于那些实时的频域分析，这些计算开销都是非常昂贵的，不过幸运的是，一个DFT的优化算法很快被开发出来，我们称之为快速傅里叶变换（Fast Fourier Transformation），当然，其本质上仍然是DFT只是算法上进行了优化使它更加适合于计算机处理，不过话说归说，不给出理论证明的广告都是瞎扯淡，那么，接下来就继续看看，DFT是如何被优化的：首先第一点，我们仍然先把离散傅里叶变换对贴上来现在，我们假设对离散序列x[n]的离散傅里叶变换写作就可以写成进一步变形，得到后面两个式子是不是非常眼熟，没错，一个DFT变换变成了2个DFT变换，不同的是，后面序列的长度只有前面的一半，为了保证后面两个DFT变换成立，k的范围也应随之变为了【0，二分之n】既然，前半部分可以变为两个DFT变换，那么后半部分呢，我们使用k+来表示离散序列后半部分，那么，DFT就变为了进一步变形得到根据欧拉公式，因为所以可以看到，仅仅是多了一个，其它都与之前的推导相同，那么这也就是为什么我们要将离散序列分为奇数列与偶数列的原因了，偶数列不变，奇数列多个负号，于是，公式变为了那么，公式总结为那么，计算的复杂度就从,只要n的值大于2，计算的时间复杂度无疑是降低了，我们进一步对多项式进行分解，直到最后仅剩下2个离散点的傅里叶变换，那么，它的复杂度将会降低至，同时，这也就要求我们的离散序列项的个数必须是2的整数次幂，因此，这种快速傅里叶变换又叫做基2快速傅里叶变换，当然同样的，也有其它基底的快速傅里叶变换，但这里就不做更多的讨论了。快速傅里叶逆变换快速傅里叶逆变换完全可按照正变换的过程进行推导，实际上就是换汤不换药，根据IDFT公式同样的，我们将离散序列X[n]分成两组，例如X[0],X[2],X[4],x[6]……x[2i]为偶数序列，记为x[2i],，而将X[1],X[3],X[5]….X[2i+1]记为X[2i+1],意思是奇(odd)序列那么，公式可以写成进一步变形，得到后式与DFT正变换，仅仅只是W指数正负的不同，我们只需要修改下符号就可以了。按照上面步骤同样推导出（因推导步骤一样，推导过程略）：使用C语言编写DFT/IDFT代码回到离散傅里叶变换对我们需要先将它们变成容易编码的格式，首先是正变换利用欧拉公式，变为然后是逆变换利用欧拉公式，变形为首先因为频域涉及到复数运算（输入信号一般为实信号）因此我们需要复数结构体complex,现定义结构体及复数的加乘运算函数之后定义DFT运算函数void DFT(complex x[],complex X[],int N);其中，x[]为输入时域信号，X[]为输出频域信号，N为时域信号的样本点数同时定义DFT逆变换函数void IDFT(complex x[],complex X[],int N);其中，X[]为输入频域信号，x[]为输出时域信号，N为频域信号的样本点数实际上通过离散傅里叶变换后频域的共轭对称性，我们仅仅需要计算前半部分就可以了，这个优化操作若读者有兴趣可以自己完成。使用C语言编写FFT/IFFT代码根据FFT的公式有那么，我们就可以采用迭代的方式，将输入序列不断地拆分重组，直到它变为2点的DFT，代码如下：void FFT_Base2(_IN _OUT complex x[],int N)
{
int exbase,exrang,i,j,k;
complex excomplex,Wnk,cx0,cx1;
if (N>>2)
{
// x[] 4 base odd/even Sort
exbase=1;
exrang=0;
while (exrang<N)
{
exrang=exbase<<2;
for (i=0;i<N/exrang;i++)//for each token
{
for (j=0;j<exbase;j++)//for each atom in token
{
excomplex=x[exrang*i+exbase+j];
x[exrang*i+exbase+j]=x[exrang*i+exbase*2+j];
x[exrang*i+exbase*2+j]=excomplex;
}
}
exbase<<=1;
}
FFT_Base2(x,N>>1);
FFT_Base2(x+(N>>1),N>>1);
for(k=0;k<N>>1;k++)
{
Wnk.re=(float)cos(-2*__PI*k/N);
Wnk.im=(float)sin(-2*__PI*k/N);
cx0=x[k];
cx1=x[k+(N>>1)];
x[k]=complexAdd(cx0,complexMult(Wnk,cx1));
Wnk.re=-Wnk.re;
Wnk.im=-Wnk.im;
x[k+(N>>1)]=complexAdd(cx0,complexMult(Wnk,cx1));
}
}
else
{
//2 dot DFT
cx0=x[0];
cx1=x[1];
x[0]=complexAdd(cx0,cx1);
cx1.im=-cx1.im;
cx1.re=-cx1.re;
x[1]=complexAdd(cx0,cx1);
}
}
void FFT(_IN complex x[],_OUT complex X[],int N)
{
int size=1;
complex *i_px;
while((size<<=1)<N);
i_px=(complex *)malloc(sizeof(complex)*size);
memset(i_px,0,sizeof(complex)*N);
memcpy(i_px,x,sizeof(complex)*N);
FFT_Base2(i_px,size);
memcpy(X,i_px,sizeof(complex)*N);
free(i_px);
}
FFT的理论则依据下面两个公式，就不再复述了：C语言代码如下：void IFFT_Base2(_IN _OUT complex X[],int N)
{
int exbase,exrang,i,j,n;
complex excomplex,Wnnk,cx0,cx1;
if (N>>2)
{
// x[] 4 base odd/even Sort
exbase=1;
exrang=0;
while (exrang<N)
{
exrang=exbase<<2;
for (i=0;i<N/exrang;i++)//for each token
{
for (j=0;j<exbase;j++)//for each atom in token
{
excomplex=X[exrang*i+exbase+j];
X[exrang*i+exbase+j]=X[exrang*i+exbase*2+j];
X[exrang*i+exbase*2+j]=excomplex;
}
}
exbase<<=1;
}
IFFT_Base2(X,N>>1);
IFFT_Base2(X+(N>>1),N>>1);
for(n=0;n<N>>1;n++)
{
Wnnk.re=(float)cos(2*__PI*n/N);
Wnnk.im=(float)sin(2*__PI*n/N);
cx0=X[n];
cx1=X[n+(N>>1)];
X[n]=complexAdd(cx0,complexMult(Wnnk,cx1));
Wnnk.re=-Wnnk.re;
Wnnk.im=-Wnnk.im;
X[n+(N>>1)]=complexAdd(cx0,complexMult(Wnnk,cx1));
}
}
else
{
//2 dot IDFT
cx0=X[0];
cx1=X[1];
X[0]=complexAdd(cx0,cx1);
cx1.im=-cx1.im;
cx1.re=-cx1.re;
X[1]=complexAdd(cx0,cx1);
}
}
void IFFT(_IN complex X[],_OUT complex x[],int N)
{
int size=1,i;
complex *i_px;
while((size<<=1)<N);
i_px=(complex *)malloc(sizeof(complex)*size);
memset(i_px,0,sizeof(complex)*N);
memcpy(i_px,X,sizeof(complex)*N);
IFFT_Base2(i_px,size);
// 1/N operate
for (i=0;i<N;i++)
{
i_px[i].re/=N;
i_px[i].im/=N;
}
memcpy(x,i_px,sizeof(complex)*N);
free(i_px);
}
信号与采样《淮南子·说山训》中有“见一叶落而知岁之将暮。”, 宋朝《文录》则引曰：“山僧不解数甲子，一叶落知天下秋。”直到今天的“一沙一世界,一花一天堂.双手握无限,刹那是永恒”都在传达着局部概括全局，一刻即是永恒的思想。如果是一个充满文艺气息的理科青年，应该很容易就能发出这样的感慨：我们都是宇宙的一份子，我们是宇宙的缩影，即便微不足道，但冥冥之中我们也是世界不可或缺的一环。但如果是一个理科大神，这会儿可要费点脑子了，大神抱着心爱的四路泰坦32G内存1T的PCIE SSD和ryzen 18x的电脑，琢磨着如何把sin（x）的波形图像完整地存入电脑，很快的大神发现即使它再把电脑容量升级一倍，他也无法存储sinx的完整波形，哪怕是一段他也办不到:首先sin(x)是周期函数，它的区间是,很遗憾，就算把整个地球的沙子都做成内存颗粒，也没有办法存储一个无穷大的量。那么存储一个周期如何，很遗憾，sin(x)是连续的函数，就算是一个周期也包含着无穷多个幅度信息，就算把整个银河系的沙子都做成内存颗粒，也无法存储无穷多的幅度信息。不过很快的，大神找到了门道，既然无法存储完整的波形，那么存储当中一些关键的点总能办得到很快，大神就找到了一个周期内的波峰和波谷,如图g.1很快，波峰和波谷的水平间距刚好是周期的一半，而波峰与x轴的垂直距离刚好就是波幅。那么最终大神用2个离散的点表示了sin(x)，并且我们也知道了，假设sin(x)的频率是n，那么我们至少要以2n的频率进行采样才能够还原出原始的波形。那么，波形函数如何以数学的形式对一个时域点进行采样呢，实际上前辈们早就定义了一个理想的函数这个函数的特点是，它仅在而x在其他的值它的定义是，的面积始终为1，因此在它的宽度无限小时，高度就变得无穷大，因此在x为0处它的值为无穷大，实际上在信号处理中，我们常常规定而非无穷大以便于我们的处理因此我们就可以得到一个采样的方程，例如对sin(x)的x=0点进行采样那么我们就可以用来表示，如果对x=2进行采样，我们只需要对进行右移位处理就可以了当然假设我们需要对多个点进行采样，那么我们完全可以写成这种形式(对sinx整数值采样)在奥本海姆的《信号与系统》7.1.1中有如下的推论设对一信号以T间隔取样，那么就有如下数学式根据卷积性质时域内的相乘对应于频域内的卷积，那么就有(P(j)为冲击串频域函数)根据公式（《信号与系统》例4.8推论,为基波频率或叫频域间距）及（*为卷积符号）于是有上式说明是频率为的周期函数，由一组移位的X(jw)叠加而成，在幅度标以的变化，当大于频带宽度的一半时（）频带不会发生堆叠，反之会发生堆叠（引用信号与系统 7.1.1 冲击串采样，详细推论请查阅书籍）。实际上上面的推论可以用一个简单的话来说：表示一个持续期为T且最高频率为W的时间函数，有2TW的样本个数就足够了。实际上这句话说的内容就是著名的香农采样定律，看！那些物理学家，数学家和诗人都是虽然表达方式不一样，但他们表达的东西常常都是一个意思，不管怎么说，他们都是会玩的。钢琴与按键识别录音频率与采样“未见其人先闻其声”说的就是靠声音识别人的道理，从声学的角度来说，各个人发出的嗓音也是各有不同的，我们大脑进化出了自动筛分频率的功能，因此尽管我们见不到真人，仍然可以从他说话的声音分辨出他，在一个嘈杂的环境中，我们也可以在分辨出那个人在说什么，看来，人脑也是一个实时的滤波系统。当然，人并不是能听到所有的声音，声音归根结底仍然是波在介质中传播，人只能听到20~20000HZ的声波，因此，低于20HZ的波我们叫次声波，高于20000HZ的波叫超声波，根据上一章节所属的香农采样定律，如果我们要记录一段声波例如演唱会，那么我们的采样频率至少就应该是40000HZ实际上大部分的MP3，WAV，OGG等音乐媒体文件，使用的采样率是44100HZ，如果一个样本点用16bits也就是2字节来计算的话，每分钟大约就需要5M字节的数据量。在多媒体的音乐文件中，我们一般需要以下几个参数，来确定多媒体文件的声音是如何采样并如何播放的：1. SampleRate 采样速率2. Channel 声道数量3. BitPerSample 每个样本的位数依靠这几个参数，我们就可以播放音乐文件了PCM码流与WAV文件格式PCM 脉冲编码调制是Pulse Code Modulation的缩写，简单来说就是抽样、量化和编码，也就是上章节分别对应的采样速率、每个样本的位数（量化）的概括了，它相当于信号的RAW（原始数据格式）。WAV文件（波形文件）应该是最接近这种原始格式的文件了，它没有对数据额外的压缩，由一个文件头构成后，剩下的都是原始的PCM数据，对播放参数进行配置后将它写入声卡就可以直接播放出声音来了。WAV的文件头如下（C++ 引用PainterEngine Audio代码）struct WaveHeader
{
pt_uchar riff[4];
//data exchange flag
pt_dword size;
//filesize-header
pt_uchar wave_flag[4];
//wave
pt_uchar fmt[4];
//fmt
pt_dword fmt_len;
//
pt_word tag;
//format
pt_word channels;
//
pt_dword samp_freq;
//sample rate
pt_dword byte_rate;
//samplerate*bytes of per sample
pt_word block_align;
//channles * bit_samp / 8
pt_word bit_samp;
//bits per sample
};
其中，最关键的几个参数是channels，samp_freq和bit_samp。读取信息头之后，往后搜索data字符串，之后的数据都是PCM码流了。WAV文件与频谱分析WavFreq是由C++编写，UI框架使用Qt 4.8.6，播放接口使用DirectSound的一款简易的声波频谱分析器（因此你需要安装DirectX SDK）。每隔0.1秒，它会对声波取样4096个点并做FFT后将频谱图显示。作为本文的附录程序，WavFreq的源代码你可以在附件中找到，同时源码遵循GPLv3开源协议。软件界面如图g.1然后打开一个WAV波形文件，如图g.2打开后将显示这个波形文件的一些基本信息，如图g.3图g.3点击菜单上的Start Analyse，观察波形频谱，如图g.4琴键分析与声学建模为了更好的讲解频域在声学分析上举足轻重的作用，笔者以一段钢琴按键音作为示范，范例文件你同样可以在附录的文件当中找到。钢琴中分别标注了七个不同的音，使用1,2,3,4,5,6,7,8（dou rui mi fa so la xi do）进行标注，它们的频谱分别如下。从上面的频谱图中我们可以看出，钢琴音的频率有着明显的差别，随着音调的升高，主要的频率也随之右移动。在这里，我们使用一些非常简单的检波滤波手段来区分这个时候按下的音是哪个。首先我们先过滤掉那些并不主要的频域信号，例如在上图中，我们主要频域在150~350的区间内，同时，我们取2500000的度量值作为其阈值 以避免一些噪声的干扰。这样，钢琴按键的识别过程就简单变为了：当前频幅度量大于2500000时，在150~350HZ区间内找到幅值能量最高的点，依据该点所在频率，滤出一些关键的特征频率，并依此来判断按键的类别。具体的代码与实现你仍然可以在附件中找到，运行结果如下图所示（图h.1）被泄漏的密码如果你有看过那些间谍大片也许你对下面这个镜头场景并不陌生,例如《谍影重重》中就有这样一个情节，特工想要进入目标的一个保险库，于是他录制了目标人物的声音然后伪造了他的语音密码，成功入侵了目标的保险库。另一个比较经典的镜头是语音的识别，特工录制了目标的语音，很快就有设备将语音识别成了文字和数字，假如目标在谈话中或者是电话的按键音中泄露了密码，那么他就要倒大霉了。不管情节如何或者这个到底具不具备可行性，上面两个场景都和声波的频域分析有关，那么从声音还原出密码具备多少的可行性呢。在很多的直播录制节目中，很多的主播在众多的观众面前输入自己的账号与密码，当然因为密码是*字遮盖的，观众无法直接看到密码，但是不少的主播的键盘敲起来是非常的响的（例如青轴的机械键盘），我们可以很清楚的听到敲击键盘噼噼啪啪的声音，有没有可能通过对键盘音的分析建模来还原出密码呢。笔者认为这绝对是有可能的，首先，回车键 上档键 空格键之类的键因为其模具的不同，其发出的敲击声音的频率肯定是不同的，另外键盘经过长时间的敲打，因为磨损与杂质的不同其击键音也有可能改变，最后一点是一个人长期的使用键盘，其对各个键的击键的力度也是不一样的，并且录音设备与键盘各键的距离不一样，很可能从声音的能量再做进一步的筛选。当然，也许完全识别出密码也许有困难，但是对主播击键音长期的采样分析建立模型，极大可能的缩小筛选范围是绝对可行的。当然，更好的声音采样器（采样精度，采样频率）对分析的成功率也有更多的帮助，以此看来劣质的麦克风反而还保护了主播们的密码隐私了。笔者还未做过相关的实验，读者有兴趣的话，不妨试试.FFT2二维冲击采样如果我们把冲击函数拓展到二维的上，显而易见的，它应该满足这样的公式：还记得之前提到一维的采样函数么，我们对原函数进行了采样操作，那么拓展到二维方向，其冲击采样应该写成了这种格式如果对离散的函数进行采样，那么公式就由积分变为了累加如果对特定样本点进行采样例如只需要对二维冲击函数进行移位就行了二维傅里叶变换对那么我们很容易也将傅里叶变换推到到二维的方面上来其中，是一个大小为M*N的矩阵，当然，逆变换同样有因此，二维的傅里叶变换过程的算法就能简单的概括为：首先依次对一个二维矩阵的每一行做傅里叶变换，得出变换的结果后，再对行变换结果的每一列做傅里叶变换。当然，为了方便计算机处理，其正变换伪代码如下1.设二维数组包含待变换的数据2.对二维数组的每一行做傅里叶变换，并将变换结果替换原数组。3.将二维数组的行与列元素互换。4.再对该二维数组的每一行做傅里叶变换，并将变换结果放回原数组。5. 再将二维数组的行与列元素互换，其结果即为二维傅里叶变换结果。用C语言实现二维傅里叶变换对二维傅里叶变换的C语言代码如下，当中引用了之前编写的一维傅里叶变换函数，当然有了之前的基础，二维傅里叶变换的代码简单的多，您可以在本文的附录中找到其完整的代码：void FFT_2(_IN complex x[],_OUT complex X[],int N)
{
for (int i=0;i<N;i++)
{
FFT(&x[i*N],&X[i*N],N);
}
//Matrix transpose
for (int cy=0;cy<N;cy++)
{
for (int cx=cy+1;cx<N;cx++)
{
complex _t=X[cy*N+cx];
X[cy*N+cx]=X[cx*N+cy];
X[cx*N+cy]=_t;
}
}
for(int i=0;i<N;i++)
{
FFT(&X[i*N],&X[i*N],N);
}
//Matrix transpose again
for (int cy=0;cy<N;cy++)
{
for (int cx=cy+1;cx<N;cx++)
{
complex _t=X[cy*N+cx];
X[cy*N+cx]=X[cx*N+cy];
X[cx*N+cy]=_t;
}
}
}
void IFFT_2(_IN complex X[],_OUT complex x[],int N)
{
//Matrix transpose
for (int cy=0;cy<N;cy++)
{
for (int cx=cy+1;cx<N;cx++)
{
complex _t=X[cy*N+cx];
X[cy*N+cx]=X[cx*N+cy];
X[cx*N+cy]=_t;
}
}
for(int i=0;i<N;i++)
{
IFFT(&x[i*N],&X[i*N],N);
}
//Matrix transpose again
for (int cy=0;cy<N;cy++)
{
for (int cx=cy+1;cx<N;cx++)
{
complex _t=X[cy*N+cx];
X[cy*N+cx]=X[cx*N+cy];
X[cx*N+cy]=_t;
}
}
for (int i=0;i<N;i++)
{
IFFT(&X[i*N],&X[i*N],N);
}
}
数字图像的相位与频谱频谱能量与相位相信读者从高中物理中经常看的到的能量公式，常常带有平方关系，例如，能量是质量与速度的平方关系，能量是电压的平方与电阻的关系。在图像中，我们常常也使用平方关系来表示“能量”这一关系，但这并不是说这张图确实带有物理上的多少能量，它更像一种比喻，就像古时群朝大臣跪地大喊吾皇万岁万岁万万岁一样，或者避讳称皇帝叫“万岁”，但至少到今天，也没有哪个皇帝有那能耐活得到一万岁的，这顶多只是口头上的一个溜屁拍马，或者说叫起来方便。图像上的能量也是类似这种关系，它更像是某种量化，并不是指它具体带有的能量，不然得话，如果你想炸掉某样东西，你只要往他们那寄一张白纸就可以了（毕竟白色所带的能量最高），更恐怖的是，要是你敢撕掉你的作业本，它甚至会引起一场爆炸。回到正题，在图像的频谱中，我们应该如何表示这种能量的关系呢，显然的，傅里叶变换的结果是一个复信号，如果求其能量，仅需要对其实部与虚部求平方相加，当然，在频谱图中我们更多的是取复数的模。也就是对其能量进行一次开平方。当然，仅仅有频谱是无法表示一个复信号的，因此我们还需要引入相位的概念其值是虚部除以实部的arctan值，当然，这也就表明它的范围是[]。最后仍然值得提到的一个细节是，二维傅里叶变换的结果是呈现中心共轭对称的。（你仍然可以套用之前一维傅里叶变换的证明方法来证明二维的共轭对称性），证明过程在本文就不再复述了，但由共轭对称性我们就很容易推到出频域图的中心对称性（直流分量部分除外部分）。图像频域与相位有什么物理意义么？实际上《数字图像处理》中有章节专门讨论这个问题，总结为图像频域表现的是灰度信息，而相位则是位置信息，频域的低频表示图像基本的灰度变化，例如一个图像灰度变化平缓也就意味着其频率较低，而频域的高频则体现其灰度剧烈变化的部分例如物体的边缘，星空的繁星（或者叫噪声？）。因为本文讨论的并不是如何对图像进行哪些模糊锐化…之类的操作，因此更多的信息我们这里不再进行更多的讨论，我们仅仅需要知道，图像的能量主要分布在了低频当中，对频域不过分的操作并不会特别影响图像的视觉效果就可以了。FFT2 DEMO程序演示有了上述的理论基础，对图像进行二维离散傅里叶变换的过程也就水到渠成了，笔者编写了DEMO程序，你可以在本文的附录中找到这个DEMO程序和它的代码，在该程序中，首先打开的界面如图所示（i.1）选择File-àOpen来打开一个图片文件，如图i.2在这里，就以笔者参照动漫《龙与虎》逢坂大河所绘制的同人图为例（非专业画师，不喜勿喷），因为是彩色图片，因此，我们分别对其的RGB（Red Green blue）分量进行分离，分别对它们做二维傅里叶变换如图i.3其中，左边第一张图是原图，第二列是原图的RBG分量，第三列是其各分量的频谱图，第四列是其相位谱。FFTShift当然，在我们观察到的频谱图中，我们更希望将频谱显示的更易于观察，在上述的频谱中，其低频分量位于左上角，且因为傅里叶变换的共轭对称性，我们同样很容易推导出二维傅里叶变换是中心对称的，假设我们将低频分量移动到原点上向外为高频，那么图像就会更加的易于观察。FFTShift的算法并不复杂，它实际上仅仅是对矩阵进行分割后（分割为4部分），对对角线的两部分进行两两交换。如图i.4例如如下矩阵经过FFTShift后就变为了例如如下4x4矩阵经过fftshift后变为经过FFTShift后，DEMO的频谱与相位图亦发生相应的改变可以看到，原本分散于四个角的低频信号经过移位到中心后，变得更加的易于观察了。对数变换有句话叫真理往往掌握在少数人的手中，在图像的二维傅里叶变换的频谱能量图中，从上个章节我们很容易观察到，频谱的能量往往在中心才较为的明显（低频域），显然的，图像的能量往往集中在低频当中。因为低频与高频的能量差距过大，为了便于频域图像的观察，我们对频域进行对数操作，并对操作的结果缩放到0-255的灰度级别中，假设表示图像频域复信号的模值（就是能量开根号），那么对数变换就是下图是经过了对数变换的频域图像，可以看到，其频域灰度变化变得更加的平缓了（图i.6）鲁棒盲水印盲水印与版权图像水印被广泛运用在了防伪，签名，标识等版权保护方面，简单来说就是防止盗图狗窃取自己的劳动成果或者将他抓个现行。毕竟网络大了什么样的人都有，被人窃取成果反而署名上其他人的名字是一件极其令人愤慨的事情。目前的水印主要是明水印（可见水印）居多，例如下图（j.1）使用的就是明水印明水印加起来简单方便并且快捷，使用Photoshop来处理应该是一件非常简单的事情，不需要太多的技术含量就可以给这张图“冠名”了当然，简单的事情往往攻击起来也很简单，稍微懂点的盗图者很容易就能把明水印删除并还原原图，实在嫌麻烦可以直接裁剪图片，很多时候裁去水印并不会对图像整体的视觉上造成多大影响（如图j.2）。当然，最烦人的是，水印常常破坏对原图的视觉享受，这在影像原画作品中，常常是不能被观众容忍的，一个作品突然冒冒失失的加上一个水印，常常导致观众看起来就像心理有个疙瘩。那么有没有水印技术，直接看不见，经过处理后就变得可见了呢。这实际就是本文另一个要讨论的技术细节，实际上这种盲水印技术在小学的动手实验上就有，最知名的应该是用米汤写字，写完后纸上是看不见的，如果要显示信息，那么就用点碘液一洗，字就变蓝色显示出来了。使用流程图来表示这一过程（如图j.3 j.4）在数字图像中，目前非常流行的一种盲水印隐写技术，是对图像像素操作的，其具体流程是，因为一个颜色具有RGB分量，一般来说，RGB三个颜色通道每个各占8位，也就是三字节，对8位的最低位做修改，对原图的影响是微乎其微的，那么，每个像素就能带有3位的信息，一张100*100的24位位图，就能够带有30000位3750字节的信息可插入，如果将水印数据插入到这里，就可以实现一个简单的盲水印了。这种盲水印技术具有一定的可行性，但是，它的缺点仍然也是非常致命的，对原图像的旋转，平移，缩放，改变色相都能非常轻易地摧毁水印。于是，具备抗攻击的盲水印技术，也就孕育而生了。基于傅里叶变换的鲁棒频域水印从之前的章节我们已经了解了如何利用二维傅里叶变换将图像变换为频域，并且我们也了解了图像的能量主要集中在低频当中，那么如果我们将水印叠加在频域的高频中，理论上对原图的视觉上并不会有过多的影响，同时，在频域中的水印散列分部在空域（也就是逆变换后的图像）的各个部分，对频域水印的破坏将会变得更加的困难，同时，频域水印对图像的裁切，旋转，平移，加噪都有一定的抗攻击能力（详细的推到可翻阅《数字图像处理》一书），因此，频域水印作为一种盲水印手段是拥有其理论依据的，对于这种具有一定抗攻击能力的水印技术，我们又称之为鲁棒水印。最后，作为讨论的细节，如何将水印叠加到频域也是应该讨论的范畴，在这里，长话短说地总结以下几点：1. 对彩色图像加水印，首先对图像的RGB颜色通道分别分离，对R，G，B三个通道的颜色分别计算频域，就和前几个章节处理的那样。2. 叠加混合的方式分为两种，称之为缩放叠加，一种为放大一种为缩小，因为二维傅里叶变换后的结果是一个复信号，因此，如果我们仅仅修改频谱能量而不影响其相位的话，应该将复信号的实部与虚部同比例放大或缩小就可以了。因此，水印也就是安装水印轮廓对对应复信号进行放大或缩小3. 峰值信噪比（PSNR），归一化相关系数（NC值）可用于判定水印对原图的影响及水印的抗干扰能力，两者都是越高越好，当然，这常常是一个相互矛盾的度量，往往抗干扰能力强也意味着对原图的干扰大。水印程序DEMOImageSigner是一款由笔者开发的图像水印程序，作为一个演示的DEMO程序，它仅仅支持256*256大小的图片，你可以在本文的附录中找到它的完整源代码，您可以通过查阅源代码，来了解频域水印的具体算法及过程，它使用C++与Qt Framework 4.8.6编写完成，遵循GPLv3开源协议。下面介绍其水印签名及各种攻击效果。1.打开程序，界面如图（k.1）2.点击Open Source Image和Open Sign Image分别加载源图与签名图片（图k.2 k.3）载入后界面如图所示（k.4）在签名控制面板，设置签名的参数（图k.5）其中，第一栏表示签名的通道，第二栏为签名水印的混合模式，一般来说，enlarge对签名图像的抗干扰能力更强，Reduce模式抗干扰模式较弱，但Enlarge模式对原图的影响较大。Power一栏表示对签名的混合能量（值越大表示水印越明显），一般来说，能量越大对原图的干扰也较大。选择完成后，点击Do Sign进行签名（图k.6）在右侧你可以看到签名后的图片，点击Save to file可以将签名后的文件保存。 最后，我们选择签名的图像，并分别使用裁剪，平移，旋转，噪声，涂抹，改变色相，缩放攻击并查看攻击后水印的保留效果。下面图k.7为原图，k.8 k.9分别为缩小及放大混合签名后的图像，可以看到视觉区别不大。同时，签名后图像的频谱如k.10 K.11（放大水印只加在蓝色通道）所示使用PhotoShop对图像进行攻击1. 裁剪攻击与其频谱k.12（蓝色通道）1. 平移攻击及其频谱1. 旋转攻击及其频谱1. 高斯噪声与频谱1. 涂抹攻击与频谱1. 改变色相攻击与其频谱1. 缩放攻击（裁剪后放大至原尺寸）与频谱后记文章到了这里也许你对当中的技术意犹未尽，也许你看的一脸云里雾里不知所云，但不管如何，《从三角函数到语音识别再到数字水印》到了这里也应该告一段落了。我一直有个愿望，把一个看上去复杂的东西，从它最原始的最简单的根基开始，抽丝剥茧，循序渐进，一步一步地深入直到它开花结果，这个过程，就犹如见证了一场破茧化蝶。因为我相信，那些深入宇宙奥妙与规律的知识，是足以让人震撼的，它并不输于任何不可估价的艺术品。遗憾的是，如今国内不论是学术还是技术都有一股浮躁之风，有人将其当做宣扬炫耀的筹码，有人甚至剽窃他人的成果与劳动，作为自己谋利的肮脏台阶，在我看来，他们都不是真正的“为学”者，我希望有更多的人，能够沉下耐心，将前辈的知识打碎发扬，将那些看起来遥不可及的智慧，分解成一步一阶的台阶，让更多的人更快更简单地得到智慧的福泽。而不是光顾着告诉大家我掌握了这个技术有多么“了不起”。我并不知道读者们阅读本文时，是否能收到我所传达的意思与愿望，本文行文时间将近三个月，当中每一个文字每一行标点每一句代码包括每一张例示图片都出自本人之手都凝聚着我的心血，但我并非圣贤，文章并不能做到句句严谨，行行精辟，相信有不少的缺漏错误还待读者们斧正，但知识理当能够有所发扬，如果您在我的字里行间最终有所感悟并了解了“她”是怎么一回事，那么，您的一句“写的不错”将是我最大的欣慰与鼓舞。最后我还想再说些什么，却有些词穷了，我将ImageSigner做更进一步的改善，让它能够支持到2次幂高宽的图像数字签名。当然，你仍然可以在附件中找到它的完整源代码与Release版本并将它运用在实际的需求中。附件下载地址：https://pan.baidu.com/s/1pL6tgjp如果您有什么相关的疑问或讨论，可以来群：689194365找我，如果想进入黑客的世界，我推荐你们观看这个专题：你想了解的炫酷白帽黑客技能都在这！行文至此，那么，就这样吧。}