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、傅里葉級數的提出及理論依據

　　在早先，拉格朗日等數學家發現某些周期函數可以由三角函數的和來表示，比如下圖中，黑色的斜線就是周期為的函數，而紅色的曲線是三角函數之和，可以看出兩者確實近似：

 基於此，傅里葉提出了任意周期函數都可以寫成三角函數之和的猜測，其公式提出思想基本如下：

• 首先，根據周期函數的定義，常數函數是周期函數，周期為任意實數。所以，分解里面得有一個常數項。

• 任意函數可以分解和奇偶函數之和，$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]，為了適應函數的值域，正余弦級數必然會有各自不同的系數a_n和b_n

至此，我們就得到了傅里葉級數的猜想式，即

$$
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,a_n，b_n和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 &center_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;
}

　　結果如下：從左到右依次是原圖、傅里葉譜、功率譜、相位譜、相位譜傅里葉反變換結果、功率譜傅里葉反變換結果、傅里葉譜反變換結果

 

 

 

 

 

 

 

 　　

通過結果可以發現，相位譜決定着圖像結構的位置信息，功率譜決定着圖像結構的整體灰度分布的特性，如明暗、灰度變化趨勢等，反應的是圖像整體上各個方向上的頻率分量的相對強度。

 

 

 

 

參考資料：

下面這個方波分解示意圖用matlab怎么畫？

如何理解傅里葉變換公式？

傅里葉頻域，復數域，沖激函數，香農采樣（不介紹公式-只介紹是啥）另一種思維

DFT

opencv學習（十五）之圖像傅里葉變換dft

[數字圖像處理]頻域濾波(1)--基礎與低通濾波器

CS425 Lab: Frequency Domain Processing

傅里葉分析與圖像處理

二維傅里葉變換是怎么進行的？

c語言數字圖像處理（六）：二維離散傅里葉變換

【數字圖像處理】傅里葉變換在圖像處理中的應用

從菜鳥到完全學會二維傅立葉在圖像處理算法中的應用【老司機教你學傅立葉】

快速傅里葉變換FFT的C程序代碼實現

2DFT and 3DFT

EE261 - The Fourier Transform and its Applications

OpenCV C++使用傅里葉譜和相角重建圖像