一、时域和频域
1.1 时域
- 时域(时间域)——自变量是时间,即横轴是时间,纵轴是信号的变化。其动态信号x(t) 是描述信号在不同时刻取值的函数。
图1是正弦波的时域图,示出了振幅与时间的关系。
在时域图中,横轴是时间,纵轴是振幅。时域图显示振幅随时间的变化,可以看出峰值振幅为5V,可以算出频率f=6 Hz。
1.2 频域
- 频域(频率域)——自变量是频率,即横轴是频率,纵轴是该频率信号的幅度,也就是通常说的频谱图。
下面是图文讲解:
图2是图1中正弦波的频域图。
在频域图中,横轴是频率,纵轴是峰值振幅。频域图仅仅示出峰值振幅与频率,而不显示振幅随时间的变化。
从频域图可以看出,正弦波的频率为6Hz,这个6Hz的正弦波的峰值振幅为5V 。
频域图的优点是,从频域图中,可以一眼看出正弦波的频率和峰值振幅,整个正弦波在频域图上只是一个立柱,立柱的位置显示了正弦波的频率,立柱的高度显示了正弦波的峰值振幅
贯穿时域和频域的方法之一,就是傅立叶分析,傅立叶分析又分为两个部分:傅立叶级数和傅立叶变换。
二、傅立叶级数
傅立叶级数指出任何周期函数都可以看作不同振幅,不同相位的正弦波的叠加。
对比傅立叶变换:傅立叶变换指出非周期的函数(函数曲线下的面积是有限的)也可以用正弦或余弦乘以加权函数的积分来表示。
说的过程大概是这样子的:
在傅立叶级数中要介绍两个概念:频谱(幅度谱),相位谱。
有了这两个东西之后我们就可以更容易理解把周期函数拆分为各个正弦函数叠加的过程了。
2.1 频谱(幅度谱)
之前我们提到了时域和频域,从不同的“域”来看可能会产生很不一样的效果,到底有多不一样呢?先上个图来看一下。
可以看出,从时域来看,我们会看到一个近似为矩形的波,而我们知道这个矩形的波可以被差分为一些正弦波的叠加。
而从频域方向来看,我们就看到了每一个正弦波的幅值,可以发现,在频谱中,偶数项的振幅都是0,也就对应了图中的彩色直线。振幅为 0 的正弦波。
上图也动态展示了频域图像,应该可以加深理解。
2.2 相位谱
频谱只代表了一个正弦函数的幅值,而要准确描述一个正弦函数,我们不仅需要幅值,还需要相位,不同相位决定了波的位置,所以对于频域分析,仅仅有频谱(振幅谱)是不够的,我们还需要一个相位谱。
那相位谱在哪呢?
先上个图
投影点我们用粉色点来表示,红色的点表示离正弦函数频率轴最近的一个峰值,而相位差就是粉色点和红色点水平距离除以周期
将相位差画到一个坐标轴上就形成了相位谱,如上图所示。
三、傅里叶变换
傅里叶级数,在时域是一个周期且连续的函数,而在频域是一个非周期离散的函数。
傅里叶变换,则是将一个时域非周期的连续信号,转换为一个在频域非周期的连续信号。
让我们先从复数说起,下面是一个复数的一个表达形式,可以看出乘以i的效果是将数值逆时针旋转了90度
数轴与虚数轴共同构成了一个复数的平面,也称复平面。这样我们就了解到,乘虚数i的一个功能——旋转。
欧拉公式将正弦波统一成了简单的指数形式,我们来看看图像上的涵义:
欧拉公式所描绘的,是一个随着时间变化,在复平面上做圆周运动的点,随着时间的改变,在时间轴上就成了一条螺旋线。如果只看它的实数部分,也就是螺旋线在左侧的投影,就是一个最基础的余弦函数。而右侧的投影则是一个正弦函数。
欧拉公式告诉我们:正弦波的叠加,也可以理解为螺旋线的叠加在实数空间的投影。
根据欧拉公式里面的复平面,我们可以得到单个矩形波形成的螺旋图,如果你认真去看,海螺图上的每一条螺旋线都是可以清楚的看到的,每一条螺旋线都有着不同的振幅(旋转半径),频率(旋转周期)以及相位。而将所有螺旋线连成平面,就是这幅海螺图了。
上图展示了将海螺图投影到实数空间就形成了傅立叶变换的连续非周期的连续的曲线,此曲线在时域上就表现为一个矩形波的形式。
2.1 变换目的
傅里叶变换是一种信号分析方法,让我们对信号的构成和特点进行深入的、定量的研究。把信号通过频谱的方式(包括幅值谱、相位谱和功率谱)进行准确的、定量的描述。这就是傅里叶变换的主要目的。
- eg.1:
如下图所示:上面我们能看到的仅仅是一个类似正弦波的波形,其幅值在按照一定的规律变化。如何记载这个波形的信息呢?尤其是量化的记载!是很困难的。那么这个时候引入傅里叶变换就可以得到一个频谱(幅值谱),主要包括3、5、7、9次谐波,一目了然!
- eg.2:
如果地震波可被分解,找出不同的振幅和速度,那么我们可以针对地震的特定振幅和速度设计对应的抗震建筑物。
如果声波可被分解成低音和高音,我们就可以放大我们关心的部分,缩小我们不关心的部分。
如果是收音机的无线电波,那么我们就可以收听到特定频率的广播。
如果计算机数据可以用震荡波形表示,且其中包含可忽略的数据,那么就可以用傅里叶变换滤去不重要的数据。这在数据科学中被叫为“数据滤波器”。
四、 其他性质
五、OpenCV中离散傅里叶变换
离散傅里叶变换(Discrete Fourier Transform,缩写为DFT),是指傅里叶变换在时域和频域上都呈现离散的形式,将时域信号的采样变换为在离散时间傅里叶变换(DTFT)频域的采样。在形式上,变换两端(时域和频域上)的序列是有限长的,而实际上这两组序列都应当被认为是离散周期信号的主值序列。即使对有限长的离散信号做DFT,也应当对其经过周期延拓成为周期信号再进行变换。在实际应用中,通常采用快速傅里叶变换来高效计算DFT.
5.1 原理
简单来说,对一张图像使用傅里叶变换就是将它分解成正弦和余弦两部分,也就是将图像从空间域(spatial domain)转换到频域(frequency domain)。
这一转换的理论基础为:任一函数都可以表示成无数个正弦和余弦函数的和的形式。傅里叶变换就是一个用来将函数分解的工具。
二维图像的傅里叶变换可以用以下数学公式表达。
式中厂是空间域(spatial domain)值, F是频域(frequency domain)值。转换之后的频域值是复数,因此,显示傅里叶变换之后的结果需要使用实数图像(real image)加虚数图像(complex image),或者幅度图像( magitude image)加相位图像(phase image)的形式。在实际的图像处理过程中,仅仅使用了幅度图像,因为幅度图像包含了原图像的几乎所有我们需要的几何信息。然而,如果想通过修改幅度图像或者相位图像的方法来间接修改原空间图像,需要使用逆傅里叶变换得到修改后的空间图像,这样就必须同时保留幅度图像和相位图像了。
在此示例中,我们将展示如何计算以及显示傅里叶变换后的幅度图像。由于数字图像的离散性,像素值的取值范围也是有限的。比如在一张灰度图像中,像素灰度值一般在0到255之间。因此,我们这里讨论的也仅仅是离散傅里叶变换(DFT),如果需要得到图像中的几何结构信息,那你就要用到它了。
在频域里面,对于一幅图像,高频部分代表了图像的细节、纹理信息;低频部分代表了图像的轮廓信息。 如果对一幅精细的图像使用低通滤波器,那么滤波后的结果就只剩下轮廓了。这与信号处理的基本思想是相通的。如果图像受到的噪声恰好位于某个特定的“频率”范围内,则可以通过滤波器来恢复原来的图像。傅里叶变换在图像处理中可以做到图像增强与图像去噪、图像分割之边缘检测、图像特征提取、图像压缩等。
5.2 dft()函数详解
dft函数的作用是对一维或二维浮点数数组进行正向或反向离散傅里叶变换。
C++: void dft (InputArray src, OutputArray dst, int flags=0, intnonzeroRows=0)
- 第一个参数, InputArray类型的src。输入矩阵,可以为实数或者虚数。
- 第二个参数, OutputArray类型的dst。函数调用后的运算结果存在这里,其尺寸和类型取决于标识符,也就是第三个参数flags
- 第三个参数, int类型的flags。转换的标识符,有默认值0,取值可以为下表中标识符的结合。
- 第四个参数,int类型的nonzeroRows,有默认值0。当此参数设为非零时(最好是取值为想要处理的那一行的值,比如C.rows),函数会假设只有输入矩阵的第一个非零行包含非零元素(没有设置DFTINVERSE标识符),或只有输出矩阵的第一个非零行包含非零元素(置了DFT INVERSE标识符),这样的话,函数就可对其他行进行更高效的处理,以节省时间开销。这项技术尤其是在采用DFT计算矩阵卷积时非常有效。
5.3 getOptimalDFTSize()函数详解
getOptimalDFTSize函数返回给定向量尺寸的傅里叶最优尺寸大小。为了提高离散傅里叶变换的运行速度,需要扩充图像,而具体扩充多少,就由这个函数来计算得到。
C++: int getOptimalDFTSize (int vecsize)
此函数的唯一一个参数为int类型的vecsize,向量尺寸,即图片的rows,cols.
5.4 copyMakeBorder()函数详解
copyMakeBorder函数的作用是扩充图像边界。
C++: void copyMakeBorder (InputArray src, OutputArray dst, int top, int bottom, int left, int right, int borderType, const Scalar& value-Scalar() )
- 第一个参数, InputArray类型的src,输入图像,即源图像,填Mat类的对象即可。
- 第二个参数, OutputArray类型的dst,函数调用后的运算结果存在这里,即这个参数用于存放函数调用后的输出结果,需和源图片有一样的尺寸和类型,且size应该为Size(src.cols+left+right, src.rows+top+bottom )
- 接下来的4个参数分别为int类型的top, bottom, left, right,分别表示在源图像的四个方向上扩充多少像素,例如top=2, bottom=2, left=2,right=2就意味着在源图像的上下左右各扩充两个像素宽度的边界。
- 第七个参数, borderType类型的,边界类型,常见取值为BORDER CONSTANT,可参考borderlnterpolate()得到更多的细节。
- 第八个参数, const Scalar&类型的value,有默认值Scalaro),可以理解为默认值为0。当borderType取值为BORDER-CONSTANT时,这个参数表示边界值。
5.5 magnitude()函数详解
函数magnitude函数用于计算二维矢量的幅值。
C++: void magnitude (InputArray x, InputArray y, OutputArray magnitude)
- 第一个参数, InputArray类型的x,表示矢量的浮点型x坐标值,也就是实部。
- 第二个参数, InputArray类型的y,表示矢量的浮点型Y坐标值,也就是虚部。
- 第三次参数, OutputArray类型的magnitude,输出的幅值,它和第一个参数x有着同样的尺寸和类型。
下式可以表示magnitude()函数的原理:
5.6 log()函数详解
log()函数的功能是计算每个数组元素绝对值的自然对数。
C++: void log (InputArray src, OutputArray dst)
第一个参数为输入图像,第二个参数为得到的对数值。其原理如下。
5.7 normalize()函数详解
normalize的作用是进行矩阵归一化。
C++: void normalize (InputArray src, OutputArray dst, double alpha=1, double beta-0, int norm type=NORM 12, int dtype=-1, InputArraymask=noArray())
第一个参数, InputArray类型的src。输入图像,即源图像,填Mat类的对象即可。
第二个参数, OutputArray类型的dst,函数调用后的运算结果存在这里,和源图片有一样的尺寸和类型。
第三个参数, double类型的alpha,归一化后的最大值,有默认值1.
第四个参数, double类型的beta。归一化后的最大值,有默认值0
第五个参数, int类型的norm_type。归一化类型,有NORM_INF、NORM_L1、NORM_L2和NORM_MINMAX等参数可选,有默认值NORM_L2.
第六个参数, int类型的dtype,有默认值-1,当此参数取负值时,输出矩阵和src有同样的类型,否则,它和src有同样的通道数,且此时图像深度为CV_MA_DEPTH (dtype)。
第七个参数, InputArray类型的mask,可选的操作掩膜,有默认值noArray()。
5.8 示例
#include "opencv2/core/core.hpp"
#include "opencv2/imgproc/imgproc.hpp"
#include "opencv2/highgui/highgui.hpp"
#include <iostream>
using namespace cv;
int main()
{
//【1】以灰度模式读取原始图像并显示
Mat srcImage = imread("F:\\CV\\LearnCV\\files\\Maliao.jpg", 0);
if (!srcImage.data) { printf("读取图片错误,请确定目录下是否有imread函数指定图片存在~! \n"); return false; }
imshow("原始图像", srcImage);
//【2】将输入图像延扩到最佳的尺寸,边界用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));
//【3】为傅立叶变换的结果(实部和虚部)分配存储空间。
//将planes数组组合合并成一个多通道的数组complexI
Mat planes[] = { Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F) };
Mat complexI;
merge(planes, 2, complexI);
//【4】进行就地离散傅里叶变换
dft(complexI, complexI);
//【5】将复数转换为幅值,即=> 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))
magnitude(planes[0], planes[1], planes[0]);// planes[0] = magnitude
Mat magnitudeImage = planes[0];
//【6】进行对数尺度(logarithmic scale)缩放
magnitudeImage += Scalar::all(1);
log(magnitudeImage, magnitudeImage);//求自然对数
//【7】剪切和重分布幅度图象限
//若有奇数行或奇数列,进行频谱裁剪
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);
//【8】归一化,用0到1之间的浮点值将矩阵变换为可视的图像格式
//此句代码的OpenCV2版为:
//normalize(magnitudeImage, magnitudeImage, 0, 1, CV_MINMAX);
//此句代码的OpenCV3版为:
normalize(magnitudeImage, magnitudeImage, 0, 1, NORM_MINMAX);
//【9】显示效果图
imshow("频谱幅值", magnitudeImage);
waitKey();
return 0;
}