OpenCV 離散傅立葉變換

摘要： 離散傅立葉變換（Discrete Fourier Transform，縮寫為DFT），是傅立葉變換在時域和頻域上都呈離散的形式，將訊號的時域取樣變換為其DFT的頻域取樣。 對於N點序列{X[n]}(0 <= n <= N),它的離散傅立葉變換為: ...

離散傅立葉變換（Discrete Fourier Transform，縮寫為DFT），是傅立葉變換在時域和頻域上都呈離散的形式，將訊號的時域取樣變換為其DFT的頻域取樣。

對於N點序列{X[n]}(0 <= n <= N),它的離散傅立葉變換為:

dft() 函式

dft()函式的作用是對一維或二維的浮點數陣列進行正向或反向的離散傅立葉變換。

函式原型

void dft(
InoutArray src,//輸入矩陣
OutputArray dst, ///輸出矩陣
int flags = 0,//轉換的識別符號
int onozeroRows,
)
複製程式碼

第三個引數，轉換的識別符號分為:

• DFT_INVERSE 用一維或二維逆變換代替預設的正向變換。
• DFT_SCALE 縮放比例識別符號，輸出的結果都會以1 / N進行縮放。
• DFT_CMPLEX_OUTPUT、DFT_REAL_OUTPUT 進行一維或二維的陣列反變換。

返回DFT最優尺寸大小：getOptimalDFTSize()函式

void int getOptimalDFTSize(int vecsize);
複製程式碼

擴充影象邊界：copyMakeBorder()函式

void copyMakeBorder(
InputArray src,//輸入影象
OutputArray dst,//輸出影象
int top,//在影象上方擴充的畫素值
int bottom,//在影象下方擴充的畫素值
int left,//在影象左方擴充的畫素值
int right,//在影象右方擴充的畫素值
int borderType, //邊界型別·
const Scalars,
)
複製程式碼

計算二維向量的幅值：magnitude()函式

用於計算二維向量的幅值

void magnitude(
InputArray x,//表示向量的浮點型X座標值,即實部
InputArray y,//表示向量的浮點型Y座標值，即虛部
OutputArray magnitude,//輸出的幅值
)
複製程式碼

計算自然對數：log() 函式

log()函式的功能是計算每個陣列元素絕對值的自然對數

void log(
InputArray src,
OutputArray dst,
);

原理即為：

if(src(I) != 0)
log|src(I)|
else
C
複製程式碼

矩陣歸一化：normalize()函式

void normalize(
InputArray src,
OutputArray dst,
double alpha = 1,//歸一化之後的最大值，有預設值1
double beta = 0,//歸一化之後的最大值，有預設值0
int norm_type = NORM_L2,//歸一化型別
int dtype = -1,//為負數時輸出矩陣和src有同樣的型別，否則，它和src有同樣的通道數，深度為CV_MAT_DEPTH
InputArray mask=noArray(),//可選的操作掩模
)
複製程式碼

getOptimalDFTSize函式

該函式返回給定向量尺寸的傅立葉最優尺寸大小。為了提高離散傅立葉變換的執行速度，需要擴充影象。

使用dft函式計算兩個二維實矩陣卷積的示例核心片段

• 其中MulSpectrums的作用是計算兩個傅立葉頻譜的每個元素的乘法

void convolveDft(InputArray A, InputArray B, OutputArray C)
{
//初始化輸出矩陣
C.create(abs(A.rows - B.rows) + 1, abs(A.cols - B.cols) + 1, A.type);

//計算DFT變換的尺寸
dftSize.width = getOptimalDFTSize(A.cols + B.cols - 1);
dftSize.height = getOptimalDFTSize(A.rows + B.rows - 1);

//分配臨時緩衝區並初始化置0
Mat tempA(dftSize,A.type(),Scalar::all(0));
Mat tempB(dftSize,B.type(),Scalar::all(0));

//分別複製A和B到tempA和tempB的左上角
Mat roiA(tempA,Rect(0,0,A.cols,A.rows));
A.copyTo(roiA);
Mat roiB(tempB,Rect(0,0,B.cols,B.rows));
B.copyTO(roiB);

//就地操作，進行快速傅立葉變換，並將nonzeroRows 引數置為零，以進行更加快速的處理。
dft(tempA,tempA,0,A.rows);
dft(tempB,tempB,0,B.rows);

//將得到的頻譜相乘，結果存放於tempA當中
mulSpectrums(tempA,tempB,tempA);

//將結果變換為頻域，儘管結果行（result.rows）都為非零，我們只需其中的C.rows的第一行，所以採用nonzeroRows == C.rows
dft(tempA,tempA,EFT_INVERSE + EFT_SCALE,C.rows);

//將結果複製到C當中
tempA(Rect(0,0,C.cols,C.rows)).copyTo(C);

//所有的臨時快取區將被自動釋放，所以無需收尾操作
}
複製程式碼

例示程式：離散傅立葉變換

【1】載入原始影象

//【1】ui灰度模式讀取原始影象並顯示
Mat srcImage = imread("D:\\Desktop\\lena.jpg",0);
if (!srcImage.data)
{
cout << "lena圖片讀取錯誤" << endl;
return false;
}
imshow("原始影象", srcImage);
複製程式碼

【2】將影象擴大到合適的尺寸

//【2】將輸入的影象延擴到最佳的尺寸，邊界用0補充
int m = getOptimalDFTSize(I.rows);
int n = getOptimalDFTSize(I.cols);
//將新增的畫素初始化為0
Mat padded;
copyMakeBorder(I,padded,0,m - I,rows,0,n - I.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】進行離散傅立葉變換

def(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]);
Mat magnitudeImage = planes[0];
複製程式碼

【6】進行對數尺度縮放

傅立葉變換的幅度之範圍大到不合適在螢幕上顯示。高值在螢幕上顯示為白點，而低值為黑點，高低值的變化無法有效判斷。為了在螢幕上凸顯出高低變化的連續性，可以用對數尺度來代替線性尺度，公式如下：

//進行對數尺度縮放
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);
q4.copyTo(q1);
tmp.copyTo(q4);
複製程式碼

【8】歸一化

現在有了重分佈後的幅度圖，但是幅度值仍然超過可顯示範圍[0,1],我們可以使用 normalize()函式歸一化到可顯示範圍。

normalize(magnitudeImage,magnitudeImage,0,1,NORM_MINMAX);
複製程式碼

【9】顯示效果

imshow("頻譜幅值"，magnitudeImage);
複製程式碼