15、频率域滤波基础——傅里叶变换计算及应用基础
1、理解傅里叶变换
如果是理工科的学生 ，在高等数学和信号处理的课程中应该就已经学习过Fourier变换 ，但是这里还是进行一个简单的基本学习和理解，为时域转频域提供一个基础理论概念。
1、什么是傅里叶级数
周期函数的fourier级数是由正弦函数和余弦函数组成的三角级数。这里首先说结论周期为T的任意周期性函数f(t)，若满足以下迪利克雷条件：
- 在一个周期内只有有限个不连续点；
- 爱一个周期内只有有限个极大和极小值
- 积分$\int_{-\frac{-T}{2}}^{\frac{T}{2}}|f(t)| d t$存在，则f(t)可展开为如下傅里叶级数
$$
f(x)=\frac{a_{0}}{2}+\sum_{k=1}^{\infty}\left(a_{k} \cos k x+b_{k} \sin k x\right)
$$
其中
$$
\begin{aligned} a_{n} &=\frac{1}{\pi} \int_{-\pi}^{\pi} f(x) \cos n x \mathrm{d} x \quad(n=0,1,2,3, \cdots) \\ b_{n} &=\frac{1}{\pi} \int_{-\pi}^{\pi} f(x) \sin n x \mathrm{d} x \quad(n=1,2,3, \cdots) \end{aligned}
$$
式中的$w=\frac{2 \pi}{T}$称为角频率。
为了便于理解，这里先给出几个傅里叶变换拟合周期函数的示例图（关于这个图如何绘制，不是本文重点，可以查看参考资料）：
2、傅里叶级数的提出及理论依据
在早先，拉格朗日等数学家发现某些周期函数可以由三角函数的和来表示，比如下图中，黑色的斜线就是周期为%22%20aria-hidden%3D%22true%22%3E%0A%3Cg%20class%3D%22mjx-svg-mrow%22%3E%0A%3Cg%20class%3D%22mjx-svg-mn%22%3E%0A%20%3Cuse%20xlink%3Ahref%3D%22%23E1-MJMAIN-32%22%3E%3C%2Fuse%3E%0A%3C%2Fg%3E%0A%3Cg%20class%3D%22mjx-svg-mi%22%20transform%3D%22translate(500%2C0)%22%3E%0A%20%3Cuse%20xlink%3Ahref%3D%22%23E1-MJMATHI-3C0%22%3E%3C%2Fuse%3E%0A%3C%2Fg%3E%0A%3C%2Fg%3E%0A%3C%2Fg%3E%0A%3C%2Fsvg%3E)
的函数，而红色的曲线是三角函数之和，可以看出两者确实近似：
基于此，傅里叶提出了任意周期函数都可以写成三角函数之和的猜测，其公式提出思想基本如下：
- 首先，根据周期函数的定义，常数函数是周期函数，周期为任意实数。所以，分解里面得有一个常数项。
-
任意函数可以分解和奇偶函数之和，$f(x)=\frac{f(x)+f(-x)}{2}+\frac{f(x)-f(-x)}{2}=f_{\text { even }}+f_{\text { odd }}$，由正余弦函数分别为奇偶函数，且其为周期性函数，而且还具有正交的特性，所以分解里面需要由正余弦函数sin和cos
- 如下图所示，对于$\sin (n x), n \in \mathbb{N}$,虽然最小周期是$\frac{\pi}{n}$,但是其周期中都有一个周期为2$\pi$，则对于周期为T的函数，可以知道角频率$w=\frac{2 \pi}{T}$，将这些函数进行加减（这里用级数表示），就保证了得到的函数的周期也为 T。
- 又因为正余弦函数的值域只有[-1，1]，而需要表示的函数的至于变化显然不止[-1，1]，为了适应函数的值域，正余弦级数必然会有各自不同的系数an和bn
至此，我们就得到了傅里叶级数的猜想式，即
$$
f(x)=\frac{a_{0}}{2}+\sum_{n=1}^{\infty}\left(a_{n} \cos \left(\frac{2 \pi n}{T} x\right)+b_{n} \sin \left(\frac{2 \pi n}{T} x\right)\right), \frac{a_{0}}{2} \in \mathbb{R}
$$
3、傅里叶级数的系数计算
在前面我们已经提出了傅里叶级数的猜想式
$$
f(x)=\frac{a_{0}}{2}+\sum_{n=1}^{\infty}\left(a_{n} \cos \left(\frac{2 \pi n}{T} x\right)+b_{n} \sin \left(\frac{2 \pi n}{T} x\right)\right), \frac{a_{0}}{2} \in \mathbb{R}
$$
但是我们还不知道系数C,an，bn和f(x)之间的关系。
为了求解系数域函数之间的关系，这里我们进一步对假设级数进行逐步积分
- 首先，对猜想式两端在$[-\pi, \pi]$上进行积分，并在等式的右端逐项积分，利用三角函数的正交性有
$$
\int_{-\pi}^{\pi} f(x) d x=\int_{-\pi}^{\pi} \frac{a_{0}}{2} d x+\sum_{n=1}^{\infty}\left(a_{n} \int_{-\pi}^{\pi} \cos n x d x+b_{n} \int_{-\pi}^{\pi} \sin n x d x\right)=a_{0} \pi
$$
于是有：
$$
a_{0}=\frac{1}{ \pi} \int_{-\pi}^{\pi} f(x) d x
$$
- 将猜想式两端同乘cos nx,再用同样的方法求积分，利用三角函数系的正交性有：
$$
\int_{-\pi}^{\pi} f(x) \cos m x d x=\int_{-\pi}^{\pi} \frac{a_{0}}{2} \cos m x d x+\sum_{n=1}^{\infty}\left(a_{n} \int_{-\pi}^{\pi} \cos n x \cos n x d x+b_{n} \int_{-\pi}^{\pi} \cos m x \sin n x d x\right)=\int_{-\pi}^{\pi} a_{n} \cos ^{2} n x=a_{n} \pi
$$
所以我们可以得出
$$
a_{n}=\frac{1}{\pi} \int_{-\pi}^{\pi} f(x) \cos n x d x
$$
- 同理将猜想式两端同乘sin nx,再用同样的方法求积分，可以得出
$$
b_{n}=\frac{1}{\pi} \int_{-\pi}^{\pi} f(x) \sin n x d x
$$
3、傅里叶级数的指数形式
结合欧拉公式
$$
e^{i t}=\cos (t)+i \sin (t)
$$
我们有
$$
\begin{aligned} \sin \theta &=\frac{e^{i \theta}-e^{-i \theta}}{2 i} \\ \cos \theta &=\frac{e^{i \theta}+e^{-i \theta}}{2} \end{aligned}
$$
根据上式，我们可以写出傅立叶级数的另外一种形式：
$$
f(x)=\sum_{n=-\infty}^{\infty} c_{n} \cdot e^{i \frac{2 \pi n x}{T}}
$$
其中
$$
c_{n}=\frac{1}{T} \int_{x_{0}}^{x_{0}+T} f(x) \cdot e^{-i \frac{2 \pi n x}{T}} d x
$$
4、由傅里叶级数拓展到傅里叶变换
我们已经知道了周期性函数只要满足迪利克雷条件就可以转换成傅里叶级数。对于非周期函数，因为其周期T趋于无穷大，我们不能将其用傅里叶级数展开，而需要做某些修改。
若f(t)为非周期函数，则可视它为周期无限大，角频率趋于0的周期函数。此时，各个相邻的谐波频率之差$\Delta \omega=(n+1) \omega_{0}-n \omega_{0}=\omega_{0}$很小，谐波频率$n \omega_{0}$需用一个变量\(\omega\)代替（此时的\(\omega\)不同于之前的角频率）。这样傅里叶级数的指数形式即可改写为
$$
\begin{array}{l}{f(t)=\sum_{\omega=-\infty}^{\infty} a_{\omega} e^{-j w t} d t} \\ {a_{w}=\frac{\Delta \omega}{2 \pi} \int_{-\frac{T}{2}}^{\frac{T}{2}} f(t) e^{-j \omega t} d t}\end{array}
$$
两式合并有：
$$
f(t)=\sum_{\omega=-\infty}^{\infty}\left[\frac{\Delta \omega}{2 \pi} \int_{-\frac{T}{2}}^{\frac{T}{2}} f(t) e^{-j \omega t} d t\right] e^{j w t}=\frac{1}{2 \pi} \sum_{\omega=-\infty}^{\infty}\left[\int_{-\frac{T}{2}}^{\frac{T}{2}} f(t) e^{-j \omega t} d t\right] e^{j w t} \Delta \omega
$$
当T趋近于无穷时，$\Delta \omega$趋近于$d \omega$，求和式变为积分式，上面的公式可以写成
$$
f(t)=\frac{1}{2 \pi} \int_{-\infty}^{\infty}\left[\int_{-\infty}^{\infty} f(t) e^{-j \omega t} d t\right] e^{j \omega t} d t
$$
则此公式即为非周期函数的傅里叶积分形式之一。
若令
$$
F(\omega)=\int_{-\infty}^{\infty} f(t) e^{-j \omega t} d t
$$
$$
f(t)=\frac{1}{2 \pi} \int_{-\infty}^{\infty} F(\omega) e^{j \omega t} d t
$$
则我们称这两个公司为傅里叶变换对，F(w)称为f(t)的傅里叶变换，而f(t)称为傅里叶反变换。
5、傅里叶变换的离散化——单变量的离散傅里叶变换（DFT）
由于计算机处理的都是离散的变量，所以连续函数必须转换为离散值序列。这是用取样和量化来完成的。对于一个连续函数f(t)，我们希望以自变量t的均匀间隔取样，由信号与系统可知，我们可以用一个单位间隔的脉冲响应作为取样函数乘以f(t),即
$$
\tilde{f}(t)=f(t) s_{\Delta T}(t)=\sum_{n=-\infty}^{\infty} f(t) \sigma\left(t-n_{\Delta} T\right)
$$
用图像表示为：
由图可知，序列中的任意取样值$f_{k}$可用如下式子表示：
$$
f_{k}=\int_{-\infty}^{\infty} f(t) \delta(t-k \Delta T) \mathrm{d} t=f(k \Delta T)
$$
取样后的函数的傅里叶变换$\tilde{F}(\mu)$是
$$
\tilde{F}(\mu)=\mathfrak{J}\{\tilde{f}(t)\}=\mathcal{J}\left\{f(t) S_{\Delta T}(t)\right\}=\mathrm{F}(\mu) \otimes S(\mu)
$$
对于脉冲串$S_{\Delta T}(t)$是周期为$\Delta T$的周期函数，则由傅里叶级数定理有
$$
s_{\Delta T}(t)=\sum_{n=-\infty}^{\infty} c_{n} \mathrm{e}^{\mathrm{j} \frac{2 \pi n}{\Delta T} t}
$$
式中：
$$
c_{n}=\frac{1}{\Delta T} \int_{-\Delta T / 2}^{\Delta T / 2} s_{\Delta T}(t) \mathrm{e}^{-\mathrm{j} \frac{2 \pi n}{\Delta T} t} \mathrm{d} t
$$
由冲激串图可以看出
在区间$[-\Delta T / 2, \Delta T / 2]$的积分只包含位于原点的冲激，即
$$
c_{n}=\frac{1}{\Delta T} \int_{-\Delta T / 2}^{\Delta T / 2} \delta(t) \mathrm{e}^{-\mathrm{j} \frac{2 \pi n}{\Delta T}} \mathrm{d} t
$$
由冲激的采样特性
$$
\int_{-\infty}^{\infty} f(t) \delta(t) \mathrm{d} t=f(0)
$$
知
$$
\mathrm{c}_{n}=\frac{1}{\Delta T} e^{-j \frac{2 \pi n}{\Delta T} 0}=\frac{1}{\Delta T} e^{0}=\frac{1}{\Delta T}
$$
因此，傅里叶级数可展开为：
$$
s_{\Delta T}(t)=\frac{1}{\Delta T} \sum_{n=-\infty}^{\infty} \mathrm{e}^{j \frac{2 \pi n}{\Delta T} t}
$$
由傅里叶变换的指数形式可以知道冲激$\delta\left(t-t_{0}\right)$的傅里叶变换为$e^{-j 2 \pi \mu t_{0}}$，所以$S(\mu)$可以化简为：
$$
S(\mu)=\Im\left\{s_{\Delta T}(t)\right\}=\Im\left\{\frac{1}{\Delta T} \sum_{n=-\infty}^{\infty} \mathrm{e}^{j \frac{2 \pi n}{\Delta T} t}\right\}=\frac{1}{\Delta T} \mathfrak{S}\left\{\sum_{n=-\infty}^{\infty} \mathrm{e}^{j \frac{2 \pi n}{\Delta T} t}\right\}=\frac{1}{\Delta T} \sum_{n=-\infty}^{\infty} \delta\left(\mu-\frac{n}{\Delta T}\right)
$$
所以这里我们可以直接得到
$$
\begin{aligned} \tilde{F}(\mu) &=F(\mu) \star S(\mu)=\int_{-\infty}^{\infty} F(\tau) S(\mu-\tau) \mathrm{d} \tau \\ &=\frac{1}{\Delta T} \int_{-\infty}^{\infty} F(\tau) \sum_{n=-\infty}^{\infty} \delta\left(\mu-\tau-\frac{n}{\Delta T}\right) \mathrm{d} \tau=\frac{1}{\Delta T} \sum_{n=-\infty}^{\infty} \int_{-\infty}^{\infty} F(\tau) \delta\left(\mu-\tau-\frac{n}{\Delta T}\right) \mathrm{d} \tau \\ &=\frac{1}{\Delta T} \sum_{n=-\infty}^{\infty} F\left(\mu-\frac{n}{\Delta T}\right) \end{aligned}
$$
至此，我们就得到了关于原始函数变换的取样过的数据的变换$\tilde{F}(\mu)$，但是这里最初我们给定的还是连续函数。
接下来我们思考如果给定的是离散的信号的时候，则$\tilde{f}(t)$的傅里叶变换为：
$$
\tilde{F}(\mu)=\int_{-\infty}^{\infty} \tilde{f}(t) \mathrm{e}^{-\mathrm{j} 2 \pi \mu t} \mathrm{d} t
$$
将最开始我们给出的$\tilde{f}(t)$函数带入有
$$
\begin{aligned} \tilde{F}(\mu) &=\int_{-\infty}^{\infty} \tilde{f}(t) \mathrm{e}^{-\mathrm{j} 2 \pi \mu t} \mathrm{d} t=\int_{-\infty}^{\infty} \sum_{n=-\infty}^{\infty} f(t) \delta(t-n \Delta T) \mathrm{e}^{-\mathrm{j} 2 \pi \mu t} \mathrm{d} t=\sum_{n=-\infty}^{\infty} \int_{-\infty}^{\infty} f(t) \delta(t-n \Delta T) \mathrm{e}^{-\mathrm{j} 2 \pi \mu t} \mathrm{d} t \\ &=\sum_{n=-\infty}^{\infty} f_{n} \mathrm{e}^{-\mathrm{j} 2 \pi \mu n \Delta T} \end{aligned}
$$
· 由给定连续函数求解的结果我们知道，其傅里叶变换$\tilde{F}(\mu)$是周期为1$/ \Delta T$的无限连续函数，因此，我们需要在一个周期内表示函数$\tilde{F}(\mu)$，这也是DFT的基础
假设我们想在周期$\mu=0$和$\mu=1 / \Delta T$得到$\tilde{F}(\mu)$的M个等间距样本。即：
$$
\mu=\frac{m}{M \Delta T}, \quad m=0,1,2, \cdots, M-1
$$
将其带入$\tilde{F}(\mu)$,则有：
$$
F_{m}=\sum_{m}^{M-1} f_{n} \mathrm{e}^{-\mathrm{j} 2 \pi m n / M}, \quad m=0,1,2, \cdots, M-1
$$
这个表达式即为我们寻找的离散傅里叶变换。此时，给定一个由$f(t)$的M个样本组成的集合$\left\{f_{n}\right\}$，我们可以通过上式来得出一个与输入样本集合离散傅里叶变换相对应的M个复离散值的样本集合$\left\{F_{m}\right\}$。
反之，给定$\left\{F_{m}\right\}$，我们可以用傅里叶反变换复原样本集$\left\{f_{n}\right\}$
$$
f_{n}=\frac{1}{M} \sum_{m=0}^{M-1} F_{m} \mathrm{e}^{\mathrm{j} 2 \pi m n / M}, \quad n=0,1,2, \cdots, M-1
$$
6、二维傅里叶变换
有了之前的基础，我们可以将傅里叶变换拓展到二维空间,二维连续傅立叶变换为:
$$
F(u, v)=\int_{-\infty-\infty}^{\infty} \int_{-\infty}^{\infty} f(x, y) e^{-j 2 \pi(\mu x+v y)} d x d y
$$
反变换为:
$$
f(x, y)=\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} F(u, v) e^{j 2 \pi(\mu x+v y)} d u d v
$$
对于一个图像尺寸为M×N的 函数的离散傅里叶变换由以下等式给出：
$$
F(u, v)=\sum_{x=0}^{M-1} \sum_{y=0}^{N-1} f(x, y) e^{-\mathrm{j} 2 \pi(u x / M+v y / N)}
$$
给出了变换函数，我们可以得到傅里叶反变换公式为:
$$
f(x, y)=\frac{1}{M N} \sum_{u=0}^{M-1} \sum_{v=0}^{N-1} F(u, v) \mathrm{e}^{\mathrm{j} 2 \pi(u x / M+v y / N)}
$$
因为离散傅里叶变换通常是复函数，因此可用极坐标形式来表示
$$
F(u, v)=|F(u, v)| \mathrm{e}^{\mathrm{j} \phi(u, v)}
$$
其中，幅度
$$
|F(u, v)|=\left[R^{2}(u, v)+I^{2}(u, v)\right]^{1 / 2}
$$
称为傅里叶变换谱（频谱）
而
$$
\phi(u, v)=\arctan \left[\frac{I(u, v)}{R(u, v)}\right]
$$
称为相角，这里的actan必须用四象限反正切来计算。
最后，功率谱定义为：
$$
F(u, v)=|F(u, v)|^{2}=R^{2}(u, v)+I^{2}(u, v)
$$
R和I分别是F(u，v)的实部和虚部。
2、基于傅里叶变换的数字信号处理
傅里叶的意义就是给出了一个在频域解决问题的方法，本质就是解释了任意波都是由简单的正弦波组成的。
利用傅里叶变换的这种特性，我们可以将时域的信号转到频域去处理，例如一个声音谱，我们可以将其转到频率谱如下：
相信大家听歌的时候都看到过相关的动态频谱波形。这里只是将声音信号转换成了频谱给大家观察，其实在日常生活中，我们从app听歌的声音都是无噪声干扰的，如果你做过从声音采集到功放输出的声音信号处理的过程，你就会知道，在声音采集的过程中，我们噪声是难以避免的，如果你需要做人声采集，你则需要将声音的300~3.4Khz的声音频率以外的噪声滤除掉。在模电中，你可以设计带通滤波器来完成相关的滤波工作。而当我们通过AD转换直接将声音采集后，我们同样可以通过DFT变换得到声音的频谱图来进行滤波操作。这个不是我们这个系列的重点，我就不做详细讲解了。
3、图像的傅里叶变换基础
在声学领域，我们可以把傅里叶变换的结果不同频段的能量的强度。在图像领域，我们可以将傅里叶变换可以看作是数学上的棱镜，将图像基于频率分解为不同的成分。当我们考虑光时,讨论它的光谱或频率谱。然而一般来说，将一幅图像的某些特定成分与其傅里叶变换直接联系联系起来通常是不可能的，但是关于傅里叶变换的频率成分和一幅图像的空间特性间的关系还是可以做出某些一般的叙述。例如，因为频率直接关系到空间变化率，因此直观的将傅里叶变换中的频率与图像中的亮度变换模式联系起来并不困难。二且图像与数字信号不同，其傅里叶变换结果并非是一维的单一曲线，而是一个二维的平面，类比关系如下：
同理，二维的傅里叶变换还可拓展至三维：
在光学方法做傅立叶变换一书中，给出了一个相关实现实验，如图，左侧是字符3的镂空平板，将其放透镜前面，用平行光照明，透镜焦平面上将显示其二维傅立叶变换。
#include "stdafx.h"
#include <iostream>
#include <opencv2/core.hpp>
#include <opencv2/highgui.hpp>
#include <opencv2/imgproc.hpp>
using namespace std;
using namespace cv;
int main()
{
Mat I = imread("2222.jpg", IMREAD_GRAYSCALE); //读入图像灰度图
//判断图像是否加载成功
if (I.empty())
{
cout << "图像加载失败!" << endl;
return -1;
}
else
cout << "图像加载成功!" << endl << endl;
Mat padded; //以0填充输入图像矩阵
int m = getOptimalDFTSize(I.rows);
int n = getOptimalDFTSize(I.cols);
//填充输入图像I，输入矩阵为padded，上方和左方不做填充处理
copyMakeBorder(I, padded, 0, m - I.rows, 0, n - I.cols, BORDER_CONSTANT, Scalar::all(0));
Mat planes[] = { Mat_<float>(padded), Mat::zeros(padded.size(),CV_32F) };
Mat complexI;
merge(planes, 2, complexI); //将planes融合合并成一个多通道数组complexI
dft(complexI, complexI); //进行傅里叶变换
//计算幅值，转换到对数尺度(logarithmic scale)
//=> log(1 + sqrt(Re(DFT(I))^2 + Im(DFT(I))^2))
split(complexI, planes); //planes[0] = Re(DFT(I),planes[1] = Im(DFT(I))
//即planes[0]为实部,planes[1]为虚部
// 计算谱 Mag = sqrt(Re^2 + Im^2)
magnitude(planes[0], planes[1], planes[0]); //planes[0] = magnitude
Mat magI = planes[0];
magI += Scalar::all(1);
log(magI, magI); //转换到对数尺度(logarithmic scale)
//如果有奇数行或列，则对频谱进行裁剪
magI = magI(Rect(0, 0, magI.cols&-2, magI.rows&-2));
//重新排列傅里叶图像中的象限，使得原点位于图像中心
int cx = magI.cols / 2;
int cy = magI.rows / 2;
Mat q0(magI, Rect(0, 0, cx, cy)); //左上角图像划定ROI区域
Mat q1(magI, Rect(cx, 0, cx, cy)); //右上角图像
Mat q2(magI, Rect(0, cy, cx, cy)); //左下角图像
Mat q3(magI, Rect(cx, cy, cx, cy)); //右下角图像
//变换左上角和右下角象限
Mat tmp;
q0.copyTo(tmp);
q3.copyTo(q0);
tmp.copyTo(q3);
//变换右上角和左下角象限
q1.copyTo(tmp);
q2.copyTo(q1);
tmp.copyTo(q2);
//归一化处理，用0-1之间的浮点数将矩阵变换为可视的图像格式
normalize(magI, magI, 0, 1, CV_MINMAX);
imshow("输入图像", I);
imshow("频谱图", magI);
waitKey(0);
return 0;
}
程序运行结果如下（原图可以通过对上图进行截取获得）：
下面来一步步看看，各个步骤的作用。如果不进行归一化处理的结果如下图所示，整个画面发白，根本看不到任何信息。
#include "stdafx.h"
#include <iostream>
#include <opencv2/core.hpp>
#include <opencv2/highgui.hpp>
#include <opencv2/imgproc.hpp>
using namespace std;
using namespace cv;
int main()
{
Mat I = imread("2222.jpg", IMREAD_GRAYSCALE); //读入图像灰度图
//判断图像是否加载成功
if (I.empty())
{
cout << "图像加载失败!" << endl;
return -1;
}
else
cout << "图像加载成功!" << endl << endl;
Mat padded; //以0填充输入图像矩阵
int m = getOptimalDFTSize(I.rows);
int n = getOptimalDFTSize(I.cols);
//填充输入图像I，输入矩阵为padded，上方和左方不做填充处理
copyMakeBorder(I, padded, 0, m - I.rows, 0, n - I.cols, BORDER_CONSTANT, Scalar::all(0));
Mat planes[] = { Mat_<float>(padded), Mat::zeros(padded.size(),CV_32F) };
Mat complexI;
merge(planes, 2, complexI); //将planes融合合并成一个多通道数组complexI
dft(complexI, complexI); //进行傅里叶变换
//计算幅值，转换到对数尺度(logarithmic scale)
//=> log(1 + sqrt(Re(DFT(I))^2 + Im(DFT(I))^2))
split(complexI, planes); //planes[0] = Re(DFT(I),planes[1] = Im(DFT(I))
//即planes[0]为实部,planes[1]为虚部
magnitude(planes[0], planes[1], planes[0]); //planes[0] = magnitude
Mat magI = planes[0];
magI += Scalar::all(1);
log(magI, magI); //转换到对数尺度(logarithmic scale)
//如果有奇数行或列，则对频谱进行裁剪
magI = magI(Rect(0, 0, magI.cols&-2, magI.rows&-2));
//重新排列傅里叶图像中的象限，使得原点位于图像中心
int cx = magI.cols / 2;
int cy = magI.rows / 2;
Mat q0(magI, Rect(0, 0, cx, cy)); //左上角图像划定ROI区域
Mat q1(magI, Rect(cx, 0, cx, cy)); //右上角图像
Mat q2(magI, Rect(0, cy, cx, cy)); //左下角图像
Mat q3(magI, Rect(cx, cy, cx, cy)); //右下角图像
//变换左上角和右下角象限
Mat tmp;
q0.copyTo(tmp);
q3.copyTo(q0);
tmp.copyTo(q3);
//变换右上角和左下角象限
q1.copyTo(tmp);
q2.copyTo(q1);
tmp.copyTo(q2);
//归一化处理，用0-1之间的浮点数将矩阵变换为可视的图像格式
//normalize(magI, magI, 0, 1, CV_MINMAX);
imshow("输入图像", I);
imshow("频谱图", magI);
waitKey(0);
return 0;
}
如果不进行象限调整，结果如下，低频出现在图片的四个角（四角发亮），而通过调整之后将低频调整到了原点附近。
#include "stdafx.h"
#include <iostream>
#include <opencv2/core.hpp>
#include <opencv2/highgui.hpp>
#include <opencv2/imgproc.hpp>
using namespace std;
using namespace cv;
int main()
{
Mat I = imread("2222.jpg", IMREAD_GRAYSCALE); //读入图像灰度图
//判断图像是否加载成功
if (I.empty())
{
cout << "图像加载失败!" << endl;
return -1;
}
else
cout << "图像加载成功!" << endl << endl;
Mat padded; //以0填充输入图像矩阵
int m = getOptimalDFTSize(I.rows);
int n = getOptimalDFTSize(I.cols);
//填充输入图像I，输入矩阵为padded，上方和左方不做填充处理
copyMakeBorder(I, padded, 0, m - I.rows, 0, n - I.cols, BORDER_CONSTANT, Scalar::all(0));
Mat planes[] = { Mat_<float>(padded), Mat::zeros(padded.size(),CV_32F) };
Mat complexI;
merge(planes, 2, complexI); //将planes融合合并成一个多通道数组complexI
dft(complexI, complexI); //进行傅里叶变换
//计算幅值，转换到对数尺度(logarithmic scale)
//=> log(1 + sqrt(Re(DFT(I))^2 + Im(DFT(I))^2))
split(complexI, planes); //planes[0] = Re(DFT(I),planes[1] = Im(DFT(I))
//即planes[0]为实部,planes[1]为虚部
magnitude(planes[0], planes[1], planes[0]); //planes[0] = magnitude
Mat magI = planes[0];
magI += Scalar::all(1);
log(magI, magI); //转换到对数尺度(logarithmic scale)
//如果有奇数行或列，则对频谱进行裁剪
magI = magI(Rect(0, 0, magI.cols&-2, magI.rows&-2));
//重新排列傅里叶图像中的象限，使得原点位于图像中心
int cx = magI.cols / 2;
int cy = magI.rows / 2;
//Mat q0(magI, Rect(0, 0, cx, cy)); //左上角图像划定ROI区域
//Mat q1(magI, Rect(cx, 0, cx, cy)); //右上角图像
//Mat q2(magI, Rect(0, cy, cx, cy)); //左下角图像
//Mat q3(magI, Rect(cx, cy, cx, cy)); //右下角图像
////变换左上角和右下角象限
//Mat tmp;
//q0.copyTo(tmp);
//q3.copyTo(q0);
//tmp.copyTo(q3);
////变换右上角和左下角象限
//q1.copyTo(tmp);
//q2.copyTo(q1);
//tmp.copyTo(q2);
//归一化处理，用0-1之间的浮点数将矩阵变换为可视的图像格式
normalize(magI, magI, 0, 1, CV_MINMAX);
imshow("输入图像", I);
imshow("频谱图", magI);
waitKey(0);
return 0;
}
如果不进行尺度调整，可以发现对比度不高，仅仅能看到中心一个亮点，说明尺度调整后，能显示更多细节。
#include "stdafx.h"
#include <iostream>
#include <opencv2/core.hpp>
#include <opencv2/highgui.hpp>
#include <opencv2/imgproc.hpp>
using namespace std;
using namespace cv;
int main()
{
Mat I = imread("2222.jpg", IMREAD_GRAYSCALE); //读入图像灰度图
//判断图像是否加载成功
if (I.empty())
{
cout << "图像加载失败!" << endl;
return -1;
}
else
cout << "图像加载成功!" << endl << endl;
Mat padded; //以0填充输入图像矩阵
int m = getOptimalDFTSize(I.rows);
int n = getOptimalDFTSize(I.cols);
//填充输入图像I，输入矩阵为padded，上方和左方不做填充处理
copyMakeBorder(I, padded, 0, m - I.rows, 0, n - I.cols, BORDER_CONSTANT, Scalar::all(0));
Mat planes[] = { Mat_<float>(padded), Mat::zeros(padded.size(),CV_32F) };
Mat complexI;
merge(planes, 2, complexI); //将planes融合合并成一个多通道数组complexI
dft(complexI, complexI); //进行傅里叶变换
//计算幅值，转换到对数尺度(logarithmic scale)
//=> log(1 + sqrt(Re(DFT(I))^2 + Im(DFT(I))^2))
split(complexI, planes); //planes[0] = Re(DFT(I),planes[1] = Im(DFT(I))
//即planes[0]为实部,planes[1]为虚部
// 计算功率谱 Mag = sqrt(Re^2 + Im^2)
magnitude(planes[0], planes[1], planes[0]); //planes[0] = magnitude
Mat magI = planes[0];
magI += Scalar::all(1);
//log(magI, magI); //转换到对数尺度(logarithmic scale)
//如果有奇数行或列，则对频谱进行裁剪
magI = magI(Rect(0, 0, magI.cols&-2, magI.rows&-2));
//重新排列傅里叶图像中的象限，使得原点位于图像中心
int cx = magI.cols / 2;
int cy = magI.rows / 2;
Mat q0(magI, Rect(0, 0, cx, cy)); //左上角图像划定ROI区域
Mat q1(magI, Rect(cx, 0, cx, cy)); //右上角图像
Mat q2(magI, Rect(0, cy, cx, cy)); //左下角图像
Mat q3(magI, Rect(cx, cy, cx, cy)); //右下角图像
//变换左上角和右下角象限
Mat tmp;
q0.copyTo(tmp);
q3.copyTo(q0);
tmp.copyTo(q3);
//变换右上角和左下角象限
q1.copyTo(tmp);
q2.copyTo(q1);
tmp.copyTo(q2);
//归一化处理，用0-1之间的浮点数将矩阵变换为可视的图像格式
normalize(magI, magI, 0, 1, CV_MINMAX);
imshow("输入图像", I);
imshow("频谱图", magI);
waitKey(0);
return 0;
}
4、傅里叶反变换复原图像
通过傅里叶变换将图像转换到频率域之后，我们可以通过滤波器操作对图像进行噪声滤波操作，关于滤波器的操作将在下一篇文章中讲到，这里不再赘述，经过滤波操作后的图像需要通过傅里叶反变换来得到傅里叶反变换的结果，OpenCV提供了傅里叶反变换的程序，下面是一个详细示例，可以帮助你理解相关API接口及其使用方法:
#include "stdafx.h"
#include<opencv2/opencv.hpp>
#include<iostream>
using namespace std;
using namespace cv;
int main()
{
//Mat input = imread(p[1], CV_LOAD_IMAGE_GRAYSCALE);//以灰度图像的方式读入图片
//如果不知到怎么传入p[1]。可以改为
Mat input = imread("2222.jpg", 0);
imshow("input", input);//显示原图
int w = getOptimalDFTSize(input.cols);
int h = getOptimalDFTSize(input.rows);//获取最佳尺寸，快速傅立叶变换要求尺寸为2的n次方
Mat padded;
copyMakeBorder(input, padded, 0, h - input.rows, 0, w - input.cols, BORDER_CONSTANT, Scalar::all(0));//填充图像保存到padded中
Mat plane[] = { Mat_<float>(padded),Mat::zeros(padded.size(),CV_32F) };//创建通道
Mat complexIm;
merge(plane, 2, complexIm);//合并通道
dft(complexIm, complexIm);//进行傅立叶变换，结果保存在自身
split(complexIm, plane);//分离通道
// 计算谱 Mag = sqrt(Re^2 + Im^2)
magnitude(plane[0], plane[1], plane[0]);//获取幅度图像，0通道为实数通道，1为虚数，因为二维傅立叶变换结果是复数
int cx = padded.cols / 2; int cy = padded.rows / 2;//一下的操作是移动图像，左上与右下交换位置，右上与左下交换位置
Mat temp;
Mat part1(plane[0], Rect(0, 0, cx, cy));
Mat part2(plane[0], Rect(cx, 0, cx, cy));
Mat part3(plane[0], Rect(0, cy, cx, cy));
Mat part4(plane[0], Rect(cx, cy, cx, cy));
part1.copyTo(temp);
part4.copyTo(part1);
temp.copyTo(part4);
part2.copyTo(temp);
part3.copyTo(part2);
temp.copyTo(part3);
//*******************************************************************
//Mat _complexim(complexIm,Rect(padded.cols/4,padded.rows/4,padded.cols/2,padded.rows/2));
//opyMakeBorder(_complexim,_complexim,padded.rows/4,padded.rows/4,padded.cols/4,padded.cols/4,BORDER_CONSTANT,Scalar::all(0.75));
Mat _complexim;
complexIm.copyTo(_complexim);//把变换结果复制一份，进行逆变换，也就是恢复原图
Mat iDft[] = { Mat::zeros(plane[0].size(),CV_32F),Mat::zeros(plane[0].size(),CV_32F) };//创建两个通道，类型为float，大小为填充后的尺寸
idft(_complexim, _complexim);//傅立叶逆变换
split(_complexim, iDft);//结果貌似也是复数
magnitude(iDft[0], iDft[1], iDft[0]);//分离通道，主要获取0通道
normalize(iDft[0], iDft[0], 1, 0, CV_MINMAX);//归一化处理，float类型的显示范围为0-1,大于1为白色，小于0为黑色
imshow("idft", iDft[0]);//显示逆变换
//*******************************************************************
plane[0] += Scalar::all(1);//傅立叶变换后的图片不好分析，进行对数处理，结果比较好看
log(plane[0], plane[0]);
normalize(plane[0], plane[0], 1, 0, CV_MINMAX);
imshow("dft", plane[0]);
waitKey(0);
return 0;
}
5、傅里叶功率谱和相位谱对图像
我们之前都是谈到的傅里叶谱及对傅里叶谱进行操作，前面我们也讲到了，傅里叶谱可以分解为功率谱和相位谱，这里通过实验讲解一下功率谱和相位谱对应整个图像的信息。
#include "stdafx.h"
#include<opencv2/opencv.hpp>
#include<iostream>
using namespace std;
using namespace cv;
//将幅度归一，相角保持不变
void one_amplitude(Mat &complex_r, Mat &complex_i, Mat &dst)
{
Mat temp[] = { Mat::zeros(complex_r.size(),CV_32FC1), Mat::zeros(complex_i.size(),CV_32FC1) };
float realv=0.0, imaginv=0.0;
for (int i = 0; i < complex_r.rows; i++) {
for (int j = 0; j< complex_r.cols; j++) {
//cout << complex_r.at<float>(i, j) << endl;
realv = complex_r.at<float>(i, j);
imaginv = complex_i.at<float>(i, j);
float distance = sqrt(realv*realv + imaginv * imaginv);
float a = realv / distance;
float b = imaginv / distance;
temp[0].at<float>(i, j) = a;
temp[1].at<float>(i, j) = b;
}
}
merge(temp, 2, dst);//多通道合成一个通道
}
//将相角归一，幅值保持不变
void one_angel(Mat &complex_r, Mat &complex_i, Mat &dst)
{
Mat temp[] = { Mat::zeros(complex_r.size(),CV_32FC1), Mat::zeros(complex_i.size(),CV_32FC1) };
float realv = 0.0, imaginv = 0.0;
for (int i = 0; i < complex_r.rows; i++) {
for (int j = 0; j < complex_r.cols; j++) {
realv = complex_r.at<float>(i, j);
imaginv = complex_i.at<float>(i, j);
float distance = sqrt(realv*realv + imaginv * imaginv);
temp[0].at<float>(i, j) = distance / sqrt(2);
temp[1].at<float>(i, j) = distance / sqrt(2);
}
}
merge(temp, 2, dst);
}
//使用1的幅值和2的相位合并
void mixed_amplitude_with_phase(Mat &real1, Mat &imag1, Mat &real2, Mat &imag2, Mat &dst)
{
if (real1.size() != real2.size()) {
std::cerr << "image don't ==" << std::endl;
return;
}
Mat temp[] = { Mat::zeros(real1.size(),CV_32FC1), Mat::zeros(real1.size(),CV_32FC1) };
float realv1 = 0.0, imaginv1 = 0.0, realv2 = 0.0, imaginv2 = 0.0;
for (int i = 0; i < real1.rows; i++) {
for (int j = 0; j < real1.cols; j++) {
realv1 = real1.at<float>(i, j);
imaginv1 = imag1.at<float>(i, j);
realv2 = real2.at<float>(i, j);
imaginv2 = imag2.at<float>(i, j);
float distance1 = sqrt(realv1*realv1 + imaginv1 * imaginv1);
float distance2 = sqrt(realv2*realv2 + imaginv2 * imaginv2);
temp[0].at<float>(i, j) = (realv2*distance1) / distance2;
temp[1].at<float>(i, j) = (imaginv2*distance1) / distance2;
}
}
merge(temp, 2, dst);
}
//使用1的相位和2的幅值合并
void mixed_phase_with_amplitude(Mat &real1, Mat &imag1, Mat &real2, Mat &imag2, Mat &dst)
{
if (real1.size() != real2.size()) {
std::cerr << "image don't ==" << std::endl;
return;
}
Mat temp[] = { Mat::zeros(real1.size(),CV_32FC1), Mat::zeros(real1.size(),CV_32FC1) };
float realv1 = 0.0, imaginv1 = 0.0, realv2 = 0.0, imaginv2 = 0.0;
for (int i = 0; i < real1.rows; i++) {
for (int j = 0; j < real1.cols; j++) {
realv1 = real1.at<float>(i, j);
imaginv1 = imag1.at<float>(i, j);
realv2 = real2.at<float>(i, j);
imaginv2 = imag2.at<float>(i, j);
float distance1 = sqrt(realv1*realv1 + imaginv1 * imaginv1);
float distance2 = sqrt(realv2*realv2 + imaginv2 * imaginv2);
temp[0].at<float>(i, j) = (realv1*distance2) / distance1;
temp[1].at<float>(i, j) = (imaginv1*distance2) / distance1;
}
}
merge(temp, 2, dst);
}
cv::Mat fourior_inverser(Mat &_complexim)
{
Mat dst;
Mat iDft[] = { Mat::zeros(_complexim.size(),CV_32F),Mat::zeros(_complexim.size(),CV_32F) };//创建两个通道，类型为float，大小为填充后的尺寸
idft(_complexim, _complexim);//傅立叶逆变换
split(_complexim, iDft);//结果貌似也是复数
magnitude(iDft[0], iDft[1], dst);//分离通道，主要获取0通道
// dst += Scalar::all(1); // switch to logarithmic scale
// log(dst, dst);
//归一化处理，float类型的显示范围为0-255,255为白色，0为黑色
normalize(dst, dst, 0, 255, NORM_MINMAX);
dst.convertTo(dst, CV_8U);
return dst;
}
void move_to_center(Mat ¢er_img)
{
int cx = center_img.cols / 2;
int cy = center_img.rows / 2;
Mat q0(center_img, Rect(0, 0, cx, cy));
Mat q1(center_img, Rect(cx, 0, cx, cy));
Mat q2(center_img, Rect(0, cy, cx, cy));
Mat q3(center_img, Rect(cx, cy, cx, cy));
Mat tmp;
q0.copyTo(tmp);
q3.copyTo(q0);
tmp.copyTo(q3);
q1.copyTo(tmp);
q2.copyTo(q1);
tmp.copyTo(q2);
}
void fast_dft(cv::Mat &src_img, cv::Mat &real_img, cv::Mat &ima_img)
{
src_img.convertTo(src_img, CV_32FC1);
///////////////////////////////////////快速傅里叶变换/////////////////////////////////////////////////////
int oph = getOptimalDFTSize(src_img.rows);
int opw = getOptimalDFTSize(src_img.cols);
Mat padded;
copyMakeBorder(src_img, padded, 0, oph - src_img.rows, 0, opw - src_img.cols,
BORDER_CONSTANT, Scalar::all(0));
Mat temp[] = { padded, Mat::zeros(padded.size(),CV_32FC1) };
Mat complexI;
merge(temp, 2, complexI);
dft(complexI, complexI);//傅里叶变换
split(complexI, temp); //将傅里叶变换结果分为实部和虚部
temp[0].copyTo(real_img);
temp[1].copyTo(ima_img);
}
int main()
{
Mat image = imread("2222.jpg", IMREAD_GRAYSCALE);
imshow("woman_src", image);
Mat woman_real, woman_imag;
Mat sqrt_real, sqrt_imag;
fast_dft(image, woman_real, woman_imag);//woman_real实部，woman_imag虚部
fast_dft(image, sqrt_real, sqrt_imag);
//
Mat img_range, img_angle;
one_amplitude(woman_real, woman_imag, img_range);//原图像的幅值归一化
one_angel(woman_real, woman_imag, img_angle);//原图像的相位归一化
//
Mat woman_amp2sqrt_angle, sqrt_amp2woman_angle;
mixed_amplitude_with_phase(woman_real, woman_imag,
sqrt_real, sqrt_imag, woman_amp2sqrt_angle);
mixed_phase_with_amplitude(woman_real, woman_imag,
sqrt_real, sqrt_imag, sqrt_amp2woman_angle);
Mat amplitude, angle, amplitude_src;
magnitude(woman_real, woman_imag, amplitude); //计算原图像幅值
phase(woman_real, woman_imag, angle); //计算原图像相位
//cartToPolar(temp[0], temp[1],amplitude, angle);
move_to_center(amplitude);//幅值亮度集合到中心
divide(amplitude, amplitude.cols*amplitude.rows, amplitude_src);
imshow("amplitude_src", amplitude_src);
amplitude += Scalar::all(1);// switch to logarithmic scale
log(amplitude, amplitude);
normalize(amplitude, amplitude, 0, 255, NORM_MINMAX); //归一化 方便显示，和实际数据没有关系
amplitude.convertTo(amplitude, CV_8U);
imshow("amplitude", amplitude);
normalize(angle, angle, 0, 255, NORM_MINMAX); //归一化 方便显示，和实际数据没有关系
angle.convertTo(angle, CV_8U);
imshow("angle", angle);
//*******************************************************************
Mat inverse_amp = fourior_inverser(img_range);//傅里叶逆变换
Mat inverse_angle = fourior_inverser(img_angle);
Mat inverse_dst1 = fourior_inverser(woman_amp2sqrt_angle);
move_to_center(inverse_angle);
imshow("相位谱复原结果", inverse_amp);
imshow("功率谱复原结果", inverse_angle);
imshow("傅里叶谱图复原结果", inverse_dst1);
waitKey(0);
return 1;
}
结果如下：从左到右依次是原图、傅里叶谱、功率谱、相位谱、相位谱傅里叶反变换结果、功率谱傅里叶反变换结果、傅里叶谱反变换结果
通过结果可以发现，相位谱决定着图像结构的位置信息，功率谱决定着图像结构的整体灰度分布的特性，如明暗、灰度变化趋势等，反应的是图像整体上各个方向上的频率分量的相对强度。
