C++實現二維離散傅立葉變換
在上一篇文章《C++實現一維離散傅立葉變換》中,我們介紹了一維訊號傅立葉變換的公式和C++實現,並闡述了頻域幅值的意義。
一維傅立葉變換只適用於一維訊號,例如音訊資料、心腦電圖等。
在影象處理中,影象訊號具有高度和寬度兩個屬性,屬於二維空間訊號。將影象訊號從空間域轉換到頻域時,需使用二維離散傅立葉變換。因此,需要將傅立葉變換從一維推廣至二維。
二維連續傅立葉變換公式如下:
經過離散化後,上述公式變為:
,其中u=0,1,2,…,M-1,v=0,1,2,…,N-1。 (式1-1)
相應地,其逆變換為:
對式1-1稍作變換,可得:
其中,剛好是一維傅立葉變換的離散形式。故二維傅立葉變換可理解為先對影象的每一行進行一維傅立葉變換,變換的結果構成一個新的矩陣,再對新矩陣的每一列進行一維傅立葉變換,這樣經過兩次變換的結果就是二維傅立葉變換的最終結果。
也就是說,二維傅立葉變換相當於先按行變換,再按列變換(也可以先按列變換,再按行變換)。最終生成的頻域資料,就是由兩個方向的頻率強度資訊疊加而成。
我們可以用虛擬碼來實現上述邏輯。假設一維傅立葉變換函式為dft1(vec),二維訊號高度為h,寬度為w,那麼虛擬碼可寫成如下形式:for (int i=0; i<h; i++)
dft1(row[i]);
for (int i=0; i<w; i++)
dft1(column[i]);在本文中,我們暫不採用上述虛擬碼思路。我們直接按式1-1來進行程式設計實現:
其中,e-i2π(ux/M+vy/N)可轉化為三角函式cos(-2π(ux/M+vy/N))+isin(-2π(ux/M+vy/N)),i為虛數單位。轉化成三角函式形式,有利於計算機程式設計實現。只需要採用兩個for迴圈,就可求出單個F(u,v)值,外面再巢狀兩個for迴圈,即可求出所有F(u,v)值。
由於涉及複數運算,需要先定義一個複數類,並實現複數運算。
在上一篇文章《C++實現一維離散傅立葉變換》中,我們已經實現了一個複數類ComplexNumber。直接引用它即可。
(注:編譯環境是VS2013,使用MFC對話方塊模板)
二維離散傅立葉變換的標頭檔案如下:
Dft2.h由於二維傅立葉變換後得到的矩陣元素數值很大,並且包含實數和虛數。為便於觀察分析,需要將變換後的結果進行求模,然後歸一化到[0,255],以便儲存為頻譜圖。歸一化後,大部分畫素灰度較低,在頻譜圖上接近黑色,肉眼不容易察覺,因此還需要使用log函式對低灰度區域進行增強。另外,由於傅立葉變換本身具有對稱性,最終生成的頻譜圖的四個角也具有對稱性,因此,一般在完成歸一化後,還有一個把頻譜原點平移到頻譜圖中心的操作。在類CDft2的宣告中,對應的求模、歸一化和中心平移函式分別為generate_spectrum()、normalize_spectrum()和shift_spectrum_to_center()。#pragma once #include "ComplexNumber.h" #define MAX_MATRIX_SIZE 4194304 // 2048 * 2048 #define PI 3.141592653 class CDft2 { public: CDft2(void); ~CDft2(void); public: bool dft2(double IN const matrix[], int IN const width, int IN const height); // 二維離散傅立葉變換 bool idft2(LPVOID OUT *pMatrix, int OUT *width, int OUT *height); // 二維離散傅立葉逆變換 void generate_spectrum(); // 對變換結果求模,生成頻譜/幅度譜 void normalize_spectrum(); // 對頻譜進行歸一化操作 bool has_dft2_matrix(); // 是否已存有變換結果 bool is_shifted_to_center(); // 是否已將頻譜原點平移到影象中心 void clear_dft2_matrix(); // 清除已有的變換結果 void print_matrix(); // 列印變換結果 void print_spectrum(); // 列印頻譜 void shift_spectrum_to_center(); // 將頻譜原點平移到影象中心 public: CComplexNumber *m_dft2_matrix; double *m_spectrum_data; protected: bool m_has_dft_matrix; bool m_is_normalized; bool m_is_spectrum_shifted; int m_dft_matrix_height; int m_dft_matrix_width; };
實現檔案如下:
Dft2.cpp
#include "StdAfx.h"
#include "Dft2.h"
CDft2::CDft2(void)
{
m_dft2_matrix = NULL;
m_spectrum_data = NULL;
m_has_dft_matrix = false;
m_is_normalized = false;
m_is_spectrum_shifted = false;
m_dft_matrix_height = 0;
m_dft_matrix_width = 0;
}
CDft2::~CDft2(void)
{
if (m_has_dft_matrix && (NULL != m_dft2_matrix) && ((m_dft_matrix_height*m_dft_matrix_width)>0))
delete[] m_dft2_matrix;
if (NULL != m_spectrum_data)
delete[] m_spectrum_data;
}
bool CDft2::has_dft2_matrix()
{
return m_has_dft_matrix;
}
bool CDft2::is_shifted_to_center()
{
return m_is_spectrum_shifted;
}
void CDft2::clear_dft2_matrix()
{
if (m_has_dft_matrix && (NULL != m_dft2_matrix) && ((m_dft_matrix_height*m_dft_matrix_width)>0)) {
delete[] m_dft2_matrix;
m_has_dft_matrix = false;
m_dft_matrix_height = 0;
m_dft_matrix_width = 0;
}
}
void CDft2::print_matrix()
{
char msg[2560] = "11111 ";
if ((!m_has_dft_matrix) || (NULL == m_dft2_matrix) || (m_dft_matrix_height <= 0) || (m_dft_matrix_width <= 0))
return;
for (int u = 0; u<m_dft_matrix_height; u++) {
for (int v = 0; v<m_dft_matrix_width; v++) {
sprintf(msg + 6, "[%d,%d]: %lf + %lfi", u + 1, v + 1, m_dft2_matrix[v + u*m_dft_matrix_width].m_rl, m_dft2_matrix[v + u*m_dft_matrix_width].m_im);
OutputDebugStringA(msg);
}
}
}
void CDft2::print_spectrum()
{
char msg[256] = "11111 ";
if ((!m_has_dft_matrix) || (NULL == m_dft2_matrix) || (m_dft_matrix_height <= 0) || (m_dft_matrix_width <= 0) || (NULL == m_spectrum_data))
return;
for (int u = 0; u<m_dft_matrix_height; u++) {
for (int v = 0; v<m_dft_matrix_width; v++) {
sprintf(msg + 6, "[%d,%d]: %lf", u + 1, v + 1, m_spectrum_data[v + u*m_dft_matrix_width]);
OutputDebugStringA(msg);
}
}
}
// Fourier transform of 2-dimension matrix
// Param1: the input matrix to be transformed
// Param 2: width of the input matrix
// Param 3: height of the input matrix
bool CDft2::dft2(double IN const matrix[], int IN const width, int IN const height)
{
if (((width*height) <= 0) || ((width*height)>MAX_MATRIX_SIZE) || (NULL == matrix))
return false;
// to avoid memory leak, make sure to clear existing data buff before executing the transformation
clear_dft2_matrix();
m_dft2_matrix = new CComplexNumber[width*height];
CComplexNumber cplTemp(0, 0);
double fixed_factor_for_axisX = (-2 * PI) / height; // evaluate -i2π/N of -i2πux/N, and store the value for computing efficiency
double fixed_factor_for_axisY = (-2 * PI) / width; // evaluate -i2π/N of -i2πux/N, and store the value for computing efficiency
for (int u = 0; u<height; u++) {
for (int v = 0; v<width; v++) {
for (int x = 0; x<height; x++) {
for (int y = 0; y<width; y++) {
double powerX = u * x * fixed_factor_for_axisX; // evaluate -i2πux/N
double powerY = v * y * fixed_factor_for_axisY; // evaluate -i2πux/N
cplTemp.m_rl = matrix[y + x*width] * cos(powerX + powerY); // evaluate f(x) * e^(-i2πux/N), which is equal to f(x) * (cos(-i2πux/N)+sin(-i2πux/N)i) according to Euler's formula
cplTemp.m_im = matrix[y + x*width] * sin(powerX + powerY);
m_dft2_matrix[v + u*width] = m_dft2_matrix[v + u*width] + cplTemp;
}
}
}
}
// now we have the transformed vector, keep the info
m_has_dft_matrix = true;
m_dft_matrix_height = height;
m_dft_matrix_width = width;
return true;
}
// Fourier transform of 2-dimension matrix
// Param1: the input matrix to be transformed
// Param 2: width of the input matrix
// Param 3: height of the input matrix
bool CDft2::idft2(LPVOID OUT *pMatrix, int OUT *width, int OUT *height)
{
char msg[256] = "11111 ";
if ((!m_has_dft_matrix) || (NULL == m_dft2_matrix) || ((m_dft_matrix_height*m_dft_matrix_width) <= 0) || (!width) || (!height))
return false;
if (*pMatrix)
delete[] * pMatrix;
*pMatrix = (LPVOID)new double[m_dft_matrix_height*m_dft_matrix_width];
CComplexNumber cplTemp(0, 0);
CComplexNumber cplResult(0, 0);
double fixed_factor_for_axisX = (2 * PI) / m_dft_matrix_height; // evaluate i2π/N of i2πux/N, and store the value for computing efficiency
double fixed_factor_for_axisY = (2 * PI) / m_dft_matrix_width; // evaluate i2π/N of i2πux/N, and store the value for computing efficiency
for (int x = 0; x<m_dft_matrix_height; x++) {
for (int y = 0; y<m_dft_matrix_width; y++) {
for (int u = 0; u<m_dft_matrix_height; u++) {
for (int v = 0; v<m_dft_matrix_width; v++) {
double powerU = u * x * fixed_factor_for_axisX; // evaluate i2πux/N
double powerV = v * y * fixed_factor_for_axisY; // evaluate i2πux/N
cplTemp.SetValue(cos(powerU + powerV), sin(powerU + powerV));
cplResult = cplResult + m_dft2_matrix[v + u*m_dft_matrix_width] * cplTemp;
}
}
((double*)*pMatrix)[y + x*m_dft_matrix_width] = cplResult.m_rl / (m_dft_matrix_height*m_dft_matrix_width);
cplResult.SetValue(0, 0);
}
}
// now we have the inverse transformed matrix, keep the info
sprintf(msg + 6, "inverse fourier result: %lX", pMatrix);
OutputDebugStringA(msg);
*width = m_dft_matrix_width;
*height = m_dft_matrix_height;
return true;
}
void CDft2::generate_spectrum()
{
if (m_has_dft_matrix && (NULL != m_dft2_matrix) && ((m_dft_matrix_height*m_dft_matrix_width)>0)) {
if (NULL != m_spectrum_data) {
delete[] m_spectrum_data;
m_is_spectrum_shifted = false;
}
m_spectrum_data = new double[m_dft_matrix_height*m_dft_matrix_width];
for (int u = 0; u<m_dft_matrix_height; u++) {
for (int v = 0; v<m_dft_matrix_width; v++) {
double a = m_dft2_matrix[v + u*m_dft_matrix_width].m_rl * m_dft2_matrix[v + u*m_dft_matrix_width].m_rl;
double b = m_dft2_matrix[v + u*m_dft_matrix_width].m_im * m_dft2_matrix[v + u*m_dft_matrix_width].m_im;
m_spectrum_data[v + u*m_dft_matrix_width] = sqrt(a + b);
}
}
}
}
void CDft2::normalize_spectrum()
{
char msg[256] = "11111 ";
if (m_has_dft_matrix && (NULL != m_dft2_matrix) && ((m_dft_matrix_height*m_dft_matrix_width)>0) && ((NULL != m_spectrum_data))) {
// get max value
double max = 0;
for (int u = 0; u<m_dft_matrix_height; u++) {
for (int v = 0; v<m_dft_matrix_width; v++) {
if (m_spectrum_data[v + u*m_dft_matrix_width] > max) {
max = m_spectrum_data[v + u*m_dft_matrix_width];
}
}
}
// normalize
if (max <= 0) {
OutputDebugStringA("11111 Dft2D::NormalizeSpectrum() max value is incorrect! function aborts.");
return;
}
else {
for (int u = 0; u<m_dft_matrix_height; u++) {
for (int v = 0; v<m_dft_matrix_width; v++) {
m_spectrum_data[v + u*m_dft_matrix_width] = (255 * m_spectrum_data[v + u*m_dft_matrix_width]) / max;
}
}
}
m_is_normalized = true;
sprintf(msg + 6, "max: %lf", max);
OutputDebugStringA(msg);
}
else
m_is_normalized = false;
}
void CDft2::shift_spectrum_to_center()
{
char msg[256] = "11111 ";
if (m_is_spectrum_shifted)
return;
if (m_has_dft_matrix && (NULL != m_dft2_matrix) && ((m_dft_matrix_height*m_dft_matrix_width)>0) && ((NULL != m_spectrum_data))) {
double *tempData = new double[m_dft_matrix_height*m_dft_matrix_width];
memcpy(tempData, m_spectrum_data, m_dft_matrix_height*m_dft_matrix_width*sizeof(double));
//移到中心
for (int u = 0; u<m_dft_matrix_height; u++) {
for (int v = 0; v<m_dft_matrix_width; v++) {
//sprintf(msg+6, "%d, %d, tempData: %lf", u+1, v+1, tempData[v+u*m_nDftMatrixWidth]);
//OutputDebugStringA(msg);
if ((u<(m_dft_matrix_height / 2)) && (v<(m_dft_matrix_width / 2))) {
m_spectrum_data[v + u*m_dft_matrix_width] =
tempData[m_dft_matrix_width / 2 + v + (m_dft_matrix_height / 2 + u)*m_dft_matrix_width];
}
else if ((u<(m_dft_matrix_height / 2)) && (v >= (m_dft_matrix_width / 2))) {
m_spectrum_data[v + u*m_dft_matrix_width] =
tempData[(v - m_dft_matrix_width / 2) + (m_dft_matrix_height / 2 + u)*m_dft_matrix_width];
}
else if ((u >= (m_dft_matrix_height / 2)) && (v<(m_dft_matrix_width / 2))) {
m_spectrum_data[v + u*m_dft_matrix_width] =
tempData[(m_dft_matrix_width / 2 + v) + (u - m_dft_matrix_height / 2)*m_dft_matrix_width];
}
else if ((u >= (m_dft_matrix_height / 2)) && (v >= (m_dft_matrix_width / 2))) {
m_spectrum_data[v + u*m_dft_matrix_width] =
tempData[(v - m_dft_matrix_width / 2) + (u - m_dft_matrix_height / 2)*m_dft_matrix_width];
}
}
}
m_is_spectrum_shifted = true;
if (tempData)
delete[] tempData;
}
}
在CDft2實現檔案中,實現二維傅立葉變換的函式dft2()主要程式碼如下(基本是按照式1-1來實現):
m_dft2_matrix = new CComplexNumber[width*height];
CComplexNumber cplTemp(0, 0);
double fixed_factor_for_axisX = (-2 * PI) / height; // evaluate -i2π/N of -i2πux/N, and store the value for computing efficiency
double fixed_factor_for_axisY = (-2 * PI) / width; // evaluate -i2π/N of -i2πux/N, and store the value for computing efficiency
for (int u = 0; u<height; u++) {
for (int v = 0; v<width; v++) {
for (int x = 0; x<height; x++) {
for (int y = 0; y<width; y++) {
double powerX = u * x * fixed_factor_for_axisX; // evaluate -i2πux/N
double powerY = v * y * fixed_factor_for_axisY; // evaluate -i2πux/N
cplTemp.m_rl = matrix[y + x*width] * cos(powerX + powerY); // evaluate f(x) * e^(-i2πux/N), which is equal to f(x) * (cos(-i2πux/N)+sin(-i2πux/N)i) according to Euler's formula
cplTemp.m_im = matrix[y + x*width] * sin(powerX + powerY);
m_dft2_matrix[v + u*width] = m_dft2_matrix[v + u*width] + cplTemp;
}
}
}
}
逆變換函式idft2()也主要是按照前面給出的逆變換公式定義來實現。求模函式generate_spectrum()的實現邏輯:對於複數x = a+bi,其模為m = mod(x) = sqrt(a2+b2),其中sqrt代表根號。
歸一化函式normalize_spectrum()的實現邏輯:先遍歷矩陣,找到最大值max,然後對於矩陣中的每一個元素matrix[i],執行matrix[i] = 255*(matrix[i]/max),從而歸一化到[0,255]。
平移中心函式shift_spectrum_to_center()的實現邏輯:將頻譜圖平分為4個象限,將1、3象限對換,2、4象限對換。主要通過座標轉換實現。現在,我們編寫並執行一個測試執行緒,對一個包含三行四列資料的二維訊號進行傅立葉變換,以驗證上述程式碼。
DWORD WINAPI test(LPVOID lParam)
{
char msg[256] = "11111 ";
int width = 4;
int height = 3;
double mat[] = { 1, 1, 3, 2, 3, 4, 123, 154, 55, 2, 22, 233 };
// Fourier transform
CDft2 dft2;
dft2.dft2(mat, width, height);
dft2.print_matrix();
dft2.generate_spectrum();
dft2.normalize_spectrum();
dft2.print_spectrum();
// inverse Fourier transform
LPVOID pMat = NULL;
int iHeight = 0;
int iWidth = 0;
dft2.idft2(&pMat, &iWidth, &iHeight);
double *iMat = (double *)pMat;
if (((iWidth*iHeight)>0) && (NULL != iMat)) {
for (int x = 0; x<iHeight; x++) {
for (int y = 0; y<iWidth; y++) {
sprintf(msg+6, "inverse fourier result %d, %d: %lf, %lf", x+1, y+1, iMat[y+x*iWidth]);
OutputDebugStringA(msg);
}
}
delete[] iMat;
}
return 0;
}
在MFC對話方塊資源中新增一個test按鈕,在按鈕事件響應函式中新增:
::CreateThread(NULL,0, test, 0, 0, NULL);
然後編譯專案。編譯成功後,先開啟DebugView日誌觀察工具,再啟動生成的exe,點選test按鈕,可以在DebugView中看到以下日誌輸出:可以看到,逆變換的結果和原始訊號完全一致。
使用Matlab的fft2()函式對原始訊號進行變換,得到的結果也和上述變換結果一致。
因此我們的實現程式碼是有效的,輸出了正確的變換結果。
當原始訊號超過512*512時,本文給出的實現程式碼執行一次變換大約需要幾十秒,這令人難以忍受。
後續我們將介紹基於蝶形分解的快速傅立葉變換,其完成一次512*512原始訊號變換隻需要幾十毫秒。