1五、頻率域濾波基礎——傅里葉變換計算及應用基礎 c語言數字圖像處理（六）：二維離散傅里葉變換【數字圖像處理】傅里葉變換在圖像處理中的應用

時間 2019-12-14

標籤頻率濾波基礎傅里葉變換計算應用 c語言數字圖像處理二維離散數字圖像處理简体版

原文原文鏈接

一、理解傅里葉變換

　　若是是理工科的學生，在高等數學和信號處理的課程中應該就已經學習過Fourier變換，可是這裏仍是進行一個簡單的基本學習和理解，爲時域轉頻域提供一個基礎理論概念。html

一、什麼是傅里葉級數

　　周期函數的fourier級數是由正弦函數和餘弦函數組成的三角級數。這裏首先說結論週期爲T的任意週期性函數f(t)，若知足如下迪利克雷條件：ios

在一個週期內只有有限個不連續點；
愛一個週期內只有有限個極大和極小值
積分$\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}
$$app

式中的$w=\frac{2 \pi}{T}$稱爲角頻率。ide

爲了便於理解，這裏先給出幾個傅里葉變換擬合周期函數的示例圖（關於這個圖如何繪製，不是本文重點，能夠查看參考資料）：函數

二、傅里葉級數的提出及理論依據

　　在早先，拉格朗日等數學家發現某些周期函數能夠由三角函數的和來表示，好比下圖中，黑色的斜線就是週期爲 $2\pi$ 的函數，而紅色的曲線是三角函數之和，能夠看出二者確實近似：post

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

首先，根據周期函數的定義，常數函數是周期函數，週期爲任意實數。因此，分解裏面得有一個常數項。
任意函數能夠分解和奇偶函數之和，$f(x)=\frac{f(x)+f(-x)}{2}+\frac{f(x)-f(-x)}{2}=f_{\text { even }}+f_{\text { odd }}$，由正餘弦函數分別爲奇偶函數，且其爲週期性函數，並且還具備正交的特性，因此分解裏面須要由正餘弦函數sin和cosui
以下圖所示，對於$\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}
$$

三、傅里葉級數的係數計算

　在前面咱們已經提出了傅里葉級數的猜測式

$$
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
$$

三、傅里葉級數的指數形式

結合歐拉公式

$$
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
$$

四、由傅里葉級數拓展到傅里葉變換

　　咱們已經知道了週期性函數只要知足迪利克雷條件就能夠轉換成傅里葉級數。對於非周期函數，由於其週期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)稱爲傅里葉反變換。

五、傅里葉變換的離散化——單變量的離散傅里葉變換（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
$$

六、二維傅里葉變換

　　有了以前的基礎，咱們能夠將傅里葉變換拓展到二維空間,二維連續傅立葉變換爲:

$$
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)的實部和虛部。

二、基於傅里葉變換的數字信號處理

　　傅里葉的意義就是給出了一個在頻域解決問題的方法，本質就是解釋了任意波都是由簡單的正弦波組成的。

利用傅里葉變換的這種特性，咱們能夠將時域的信號轉到頻域去處理，例如一個聲音譜，咱們能夠將其轉到頻率譜以下：

　　相信你們聽歌的時候都看到過相關的動態頻譜波形。這裏只是將聲音信號轉換成了頻譜給你們觀察，其實在平常生活中，咱們從app聽歌的聲音都是無噪聲干擾的，若是你作過從聲音採集到功放輸出的聲音信號處理的過程，你就會知道，在聲音採集的過程當中，咱們噪聲是難以免的，若是你須要作人聲採集，你則須要將聲音的300~3.4Khz的聲音頻率之外的噪聲濾除掉。在模電中，你能夠設計帶通濾波器來完成相關的濾波工做。而當咱們經過AD轉換直接將聲音採集後，咱們一樣能夠經過DFT變換獲得聲音的頻譜圖來進行濾波操做。這個不是咱們這個系列的重點，我就不作詳細講解了。

三、圖像的傅里葉變換基礎

　　在聲學領域，咱們能夠把傅里葉變換的結果不一樣頻段的能量的強度。在圖像領域，咱們能夠將傅里葉變換能夠看做是數學上的棱鏡，將圖像基於頻率分解爲不一樣的成分。當咱們考慮光時,討論它的光譜或頻率譜。然而通常來講，將一幅圖像的某些特定成分與其傅里葉變換直接聯繫聯繫起來一般是不可能的，可是關於傅里葉變換的頻率成分和一幅圖像的空間特性間的關係仍是能夠作出某些通常的敘述。例如，由於頻率直接關係到空間變化率，所以直觀的將傅里葉變換中的頻率與圖像中的亮度變換模式聯繫起來並不困難。二且圖像與數字信號不一樣，其傅里葉變換結果並不是是一維的單一曲線，而是一個二維的平面，類比關係以下：

　　同理，二維的傅里葉變換還可拓展至三維：

　　在光學方法作傅立葉變換一書中，給出了一個相關實現實驗，如圖，左側是字符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;
}

四、傅里葉反變換復原圖像

　　　經過傅里葉變換將圖像轉換到頻率域以後，咱們能夠經過濾波器操做對圖像進行噪聲濾波操做，關於濾波器的操做將在下一篇文章中講到，這裏再也不贅述，通過濾波操做後的圖像須要經過傅里葉反變換來獲得傅里葉反變換的結果，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;
}

五、傅里葉功率譜和相位譜對圖像

　　咱們以前都是談到的傅里葉譜及對傅里葉譜進行操做，前面咱們也講到了，傅里葉譜能夠分解爲功率譜和相位譜，這裏經過實驗講解一下功率譜和相位譜對應整個圖像的信息。

#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;
}

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

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