我在上两篇文章「手把手教你编写傅里叶动画」、「傅里叶动画专辑欣赏」中介绍了傅里叶级数的本质以及编写了一些有趣的傅里叶动画,主要讲述了周期性函数究竟是如何一步步被分解成正余弦函数的和的。但是,不幸的是我们在工程中使用的一些函数往往会有一些非周期性函数,那么我们该如何用三角函数来描述它们呢,这就是今天我要讲述的傅里叶变换。
那么傅里叶变化在实际工程中具体有哪些应用领域呢?
傅氏变换模式识别图像压缩图像降噪我只举上面三个例子来说明,但是要明白傅里叶变换的应用范围远远不止上面几个,它在电子信息、电工学、机械工程学等领域都有重大应用,举个最简单的例子,现在手机相机中的夜景降噪算法的原理就是通过傅里叶变换来实现的。
一,从傅里叶级数到傅里叶变换
傅里叶变换与傅里叶级数唯一的区别就在于前者我们分解的是周期性函数,而后者我们研究的是非周期型函数。但是,如果我们转换一下思路,如果把一个非周期性函数的周期看成无穷大,那么我们不是可以按照傅里叶级数的思路去分解这个函数了吗?
在上文中,我们知道,傅里叶级数的公式为: 其中 是以 为周期的函数,我们令 ,那么上式化为: 再令 如果
那么再回头看 : 由于 连续变化,且 ,所以求和可以写成积分形式:
逆变换为: 其中 ,这个公式我们可以明显看出,傅里叶变换的本质就是对一个函数进行时域与频域之间的相互转换。
可是看到这里,我相信大多数人依然会觉得这到底是是什么东西,我举个简单的例子吧:桌子上放了一杯混合着糖、盐、酒等作料的液体,小明想具体知道里面的成分具体是什么,又占了多少比重,那么小明就需要能够过滤盐、糖的过滤器,而傅里叶变换就是扮演者过滤器的角色,它能够把混合液体的具体成分给分解出来。听起来不可思议吧,但不管你相不相信,傅里叶当年也正是这样考虑问题的。
例如方波函数就可由多个正余弦波函数来进行叠加得到,也许你听起来觉得不太可能,但事实就是如此,根据上篇文章,多个正弦波函数,如果我们调整适当的半径速度与旋转方向,我们就可以拟合出方波函数:
我们将上面四个正弦函数累加,我们看它们的结果:
上面的动画我只累加了四个正弦波函数,已经很接近方波函数了,而随着我们累加的越来越多,最终结果也会越来越接近我们的目标函数。但是,大多数情况下我们不需要那么多特征,只需要有限个波函数就能逼近我们的目标函数。反过来,我们也可以通过这种方法,舍弃后面不必要的特征,来近似得到目标函数。JPEG图像压缩、MP3音频技术就是这个原理,例如你用微信发一张图片,一般情况下微信会自动给这张图片进行压缩处理(除非你发送原图),这样可以在保证画面的主要特征不消失的情况下,减小流量消耗。
所以通过对方波函数的分解,其实我们知道傅里叶变换其实就是把函数通过时域转换到频域的过程。那么这玩意具体在工程中有什么应用呢?今天我主要说明一下傅里叶变换在图像处理中的应用。
二,傅里叶变换在图像处理中的应用
傅里叶变换在图像处理中有重大应用,例如图像的傅里叶降噪、JPEG图像压缩技术、模式识别等等。但是图像是二维的,二维傅里叶变换公式与一维稍有区别: 其中 是图像对应每个像素矩阵的亮度或者灰度值,这个公式实际几何意义就是将一个图像分解成若干个复平面波之和。从公式我们可以看出,二维傅里叶变换就是将图像与每个不同频率的不同方向做内积运算,也就是逐行逐列的使用一维傅里叶变换。图像经过傅里叶变换后能将其空间域转化成频域,那么这样做有什么好处呢?
研究表明,在频域图像里面,高频往往对应着图像的亮度或者灰度变化剧烈的地方,例如边缘信息、噪声等信息;而低频部分往往对应着图像中亮度变化不大的地方。通过这个原理,我们可以对图像频域中特定的频率做适当的过滤,再经过傅里叶逆变换,这样我们就能实现图像的降噪、压缩、特征提取等等。例如,如果图像的噪声刚好位于频域图中的特定频率,那么我们就能通过滤波器来对其过滤掉,这样就能实现图像的降噪操作。
1,图像的傅里叶降噪
上图是一张含有噪声的图片,我们先对其做傅里叶变换,查看一下其频域图,c++代码实现如下(请先配置一下opencv,具体方法请参考我另外一篇文章手把手教你安装OpenCV与配置环境):
//函数名称:DFT
//函数描述:对图像进行傅里叶变换,查看频域图
//返回值:void
//作者:@刘亚曦
void DFT(Mat& srcImage)
{
//将输入图像延扩到最佳的尺寸,边界用0补充
int m = getOptimalDFTSize(srcImage.rows);
int n = getOptimalDFTSize(srcImage.cols);
//将添加的像素初始化为0.
Mat padded;
copyMakeBorder(srcImage, padded, 0, m - srcImage.rows, 0, n - srcImage.cols, BORDER_CONSTANT, Scalar::all(0));
//傅立叶变换的结果(实部和虚部)分配存储空间。
Mat planes[] = { Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F) };
Mat complexI;
merge(planes, 2, complexI);
//进行就地离散傅里叶变换
dft(complexI, complexI);
//将复数转换为幅值,即=> log(1 + sqrt(Re(DFT(I))^2 + Im(DFT(I))^2))
split(complexI, planes); // 将多通道数组complexI分离成几个单通道数组,planes[0] = Re(DFT(I), planes[1] = Im(DFT(I))
//定义幅度谱和相位谱
Mat mag;
magnitude(planes[0], planes[1], mag);// planes[0] = magnitude
Mat magnitudeImage = mag;
//进行对数尺度(logarithmic scale)缩放
magnitudeImage += Scalar::all(1);
log(magnitudeImage, magnitudeImage);//求自然对数
//剪切和重分布幅度图象限
magnitudeImage = magnitudeImage(Rect(0, 0, magnitudeImage.cols & -2, magnitudeImage.rows & -2));
//重新排列傅立叶图像中的象限,使得原点位于图像中心
int cx = magnitudeImage.cols / 2;
int cy = magnitudeImage.rows / 2;
Mat q0(magnitudeImage, Rect(0, 0, cx, cy)); // ROI区域的左上
Mat q1(magnitudeImage, Rect(cx, 0, cx, cy)); // ROI区域的右上
Mat q2(magnitudeImage, Rect(0, cy, cx, cy)); // ROI区域的左下
Mat q3(magnitudeImage, Rect(cx, cy, cx, cy)); // ROI区域的右下
Mat tmp;
q0.copyTo(tmp);
q3.copyTo(q0);
tmp.copyTo(q3);
q1.copyTo(tmp);
q2.copyTo(q1);
tmp.copyTo(q2);
//归一化,用0到1之间的浮点值将矩阵变换为可视的图像格式
normalize(magnitudeImage, magnitudeImage, 0, 1, NORM_MINMAX);
imshow("频域图", magnitudeImage);
waitKey(0);
}
我们观察代码运行计算的结果,也就是频域图:
这就是上面那张图片的频域图,图中中间部分代表低频信息,周围部分代表高频信息。我们通过对其中周围的高频频谱做适当的过滤,再讲过傅里叶逆变换,这样就能消除噪声。但是高频中往往包含图像的细节,这样做的代价就是图像不可避免的出现了模糊:
2,文字识别:
在文字识别领域中,我们往往要矫正文档的方向,例如有时候我们采集的图像中的文字是倾斜的,这个时候我们就可以通过傅里叶变换来实现,我们先来看下面几行文字,也就是本文的开头一段:
按照上面的代码,我们试运行一遍:
这是在没有倾斜情况下的频域图,如果这段文字发生了倾斜:
我们再来看运行后结果的频域图:
神奇的一幕发生了,我们发现图像频域的内容(图中斜线)居然和文字的旋转几何方向相关,通过这点我们就可以计算文档的旋转角度并修正偏差。
3,模式识别:
在计算机模式识别领域中,我们往往要通过样本的特征将样本划分到一定的类别中,例如让计算机给手写字母分类,尽管我们肉眼能够很容易就判断出文字的分类,但对计算机来说绝非这么简单,但是如果转换到频域之后,计算机就能够很快识别,例如字母A:
上图是两个手写字母A经过傅里叶变换后的频域图,我们就会发现它们两个非常“像”;而对于其它两个字母,例如H与M:
明显,我们发现H与M这两个字母的频域图就没那么相似,因此对于计算机来说,频域图更能揭示一张图片的特征。
随意展示一张导图内容(所有的子节点都可以打开):