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
#pragma once
#include "ComplexNumber.h"#define MAX_MATRIX_SIZE 4194304 // 2048 * 2048
#define PI 3.141592653class 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;
};
由于二维傅里叶变换后得到的矩阵元素数值很大,并且包含实数和虚数。为便于观察分析,需要将变换后的结果进行求模,然后归一化到[0,255],以便保存为频谱图。归一化后,大部分像素灰度较低,在频谱图上接近黑色,肉眼不容易察觉,因此还需要使用log函数对低灰度区域进行增强。另外,由于傅里叶变换本身具有对称性,最终生成的频谱图的四个角也具有对称性,因此,一般在完成归一化后,还有一个把频谱原点平移到频谱图中心的操作。在类CDft2的声明中,对应的求模、归一化和中心平移函数分别为generate_spectrum()、normalize_spectrum()和shift_spectrum_to_center()。
实现文件如下:
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 transformationclear_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 efficiencydouble fixed_factor_for_axisY = (-2 * PI) / width; // evaluate -i2π/N of -i2πux/N, and store the value for computing efficiencyfor (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/Ndouble powerY = v * y * fixed_factor_for_axisY; // evaluate -i2πux/NcplTemp.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 formulacplTemp.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 infom_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 efficiencydouble fixed_factor_for_axisY = (2 * PI) / m_dft_matrix_width; // evaluate i2π/N of i2πux/N, and store the value for computing efficiencyfor (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/Ndouble powerV = v * y * fixed_factor_for_axisY; // evaluate i2πux/NcplTemp.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 infosprintf(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 valuedouble 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];}}}// normalizeif (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);}elsem_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 efficiencydouble fixed_factor_for_axisY = (-2 * PI) / width; // evaluate -i2π/N of -i2πux/N, and store the value for computing efficiencyfor (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/Ndouble powerY = v * y * fixed_factor_for_axisY; // evaluate -i2πux/NcplTemp.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 formulacplTemp.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 transformCDft2 dft2;dft2.dft2(mat, width, height);dft2.print_matrix();dft2.generate_spectrum();dft2.normalize_spectrum();dft2.print_spectrum();// inverse Fourier transformLPVOID 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原始信号变换只需要几十毫秒。
C++实现二维离散傅里叶变换相关推荐
- 第4章 Python 数字图像处理(DIP) - 频率域滤波5 - 二变量函数的傅里叶变换、图像中的混叠、二维离散傅里叶变换及其反变换
目录 二变量函数的傅里叶变换 二维冲激及其取样性质 二维连续傅里叶变换对 二维取样和二维取样定理 图像中的混叠 二维离散傅里叶变换及其反变换 二变量函数的傅里叶变换 二维冲激及其取样性质 两个连续变量 ...
- Python 二维离散傅里叶变换
Python 二维离散傅里叶变换 文章目录 Python 二维离散傅里叶变换 需要的库 计算两张图片的PSNR 二维离散傅里叶变换 二维离散傅里叶逆变换 频域平移 绘制频域图像 需要的库 import ...
- 傅里叶变换 二维离散傅里叶变换
1.介绍. DFT:(Discrete Fourier Transform)离散傅里叶变换是傅里叶变换在时域和频域上都呈离散的形式,将信号的时域采样变换为其DTFT的频域采样.在形式上,变换两端(时域 ...
- 图像傅里叶变换(二维离散傅里叶变换)
图像傅里叶变换 二维离散傅里叶变换是将图像从空间域转至频域,在图像增强.图像去噪.图像边缘检测.图像特征提取.图像压缩等等应用中都起着极其重要的作用.理论基础是任意函数都可以表示成正弦函数的线性组合的 ...
- 图像二维离散傅里叶变换、幅度谱、相位谱
clear, clc I = imread('...');F = fftshift(fft2(I)); % 对图像进行二维 DFT(fft2),并移至中心位置 magn = log(abs(F)); ...
- 二维离散傅里叶变换性质
1. 线性性质(加法定理): 2. 比例性质(相似性定理) 3. 可分离性: 4. 空间位移(位移定理): 5. 频率位移: 6. 周期性: 7. 共轭对称性: 8. 旋转不变性: 9. 平均值: 1 ...
- 在二维离散傅里叶变换中进行频谱平移(MATLAB::fft2shift)的作用
图像处理开发需求.图像处理接私活挣零花钱,请加微信/QQ 2487872782 图像处理开发资料.图像处理技术交流请加QQ群,群号 271891601 懒得自己敲文字描述了,直接摘取在一个资料上看到的 ...
- C++实现二维快速傅里叶变换(FFT)
上一篇文章里,我根据DFT公式用C++实现了二维离散傅里叶变换.但跑一张300*300的图片都要好几分钟,速度实在太慢了.我研究了下快速傅里叶变换,在网上找了一些资料,然后用C++实现了二维快速傅里叶 ...
- 快速傅里叶变换 java_二维快速傅里叶变换的java实现
图像处理与模式识别课作业,没学过信号与系统(哭晕). 恶补了两天冈萨里斯的书,写一下实现原理及过程 看了网络上很多版本的概念讲解,道理我都懂,但是在将算法迁移到傅里叶变换的实现上时,出现了一些问题.接 ...
最新文章
- 16个Linux服务器监控命令
- Waiting for Debugger
- Chapter10:观察者模式
- 新华三,定义服务器虚拟化市场新格局
- ROA(资产收益率)与TCO(总所有成本)解析
- android chrome cast,有了它 任何安卓设备瞬间变身ChromeCast
- 解决Android Studio默认AppTheme 没有lable标签,不显示等问题
- 计算机打开管理工具空白,大师解决win10系统打开设备管理器却显示空白的办法...
- xcopy和robocopy
- 候选公示!高工智能汽车金球奖第二批入围年度产品/方案亮相
- C#开发WebService实例和发布
- WinRAR 5.5 破解方法 - 自己动手, 更放心
- MySQL之redo日志
- 造一个智能语音音箱!!!太简单了【语音智能管家】
- 基于OpenCV 的车道线检测方法
- python 提取列表中长度大于3的字符串,列表中什么元素都有
- MySQL 3:MySQL数据库基本操作 DQL
- 03.规格及模板管理
- 游程编码用matlab实现代码_二值图像游程编码matlab代码
- 飞卡日常进度之山外上位机的二值化模式和灰度模式