简介

DFT的中文名为为离散傅里叶变换，二维DFT算法可以计算灰度图像的频率成分，而IDFT算法可以将频率成分还原为灰度图像。
DFT算法的公式如下：
$$\begin{gathered} F(u,v)=\frac{1}{MN}\sum _{x=0}^{M-1}\sum _{y=0}^{N-1}f(x,y)e^{-j2\pi (ux/M+vy/N)} \\ 式子中：u=0,1,2,\cdots ,M-1;v=0,1,2,\cdots ,N-1。 \end{gathered}$$
而IDFT的公式如下：
$$\begin{gathered} f(x,y)=\sum _{u=0}^{M-1}\sum _{v=0}^{N-1}F(u,v)e^{j2\pi (ux/M+vy/N)} \\ 式子中：x=0,1,2,\cdots ,M-1;y=0,1,2,\cdots ,N-1。 \end{gathered}$$

算法介绍

DFT和IDFT算法的实现

首先实现DFT算法，在C语言中没有复数的概念，所以不能直接计算原来的式子，需要借助欧拉公式$e^{jx}=\cos (x)+j\sin (x)$，将原来的式子更换为实部+虚部的形式，即一个数组存储实部，另外一个数组存储虚部。因此原式子可以转换为：
$$F(u,v)=\frac{1}{MN}\sum _{x=0}^{M-1}\sum _{y=0}^{N-1}f(x,y)[\cos 2\pi (ux/M+vy/N)-j\sin 2\pi (ux/M+vy/N)]$$
由此式子可以写出下面的代码，这个是没有经过优化的的C++代码，输入为灰度图像数组，存储实部和虚部的二级指针，以及图像的宽度和高度。因为float型数组的占用的空间非常巨大，使用数组存储的话比较浪费空间，到使用的时候再创建比较合适。并且这个函数是借助了math.h的，无法在单片机平台直接使用（下面给出了另外的方法）。

#define PI 3.1415926//圆周率的定义
//求取dft的实部和虚部
void discrete_fourier_transform(unsigned char img_gray[][Width], float **real, float **imaginary, unsigned int height, unsigned int width)
{
	double re, im, temp;
	for (unsigned int i = 0; i < height; i++)
	{
		for (unsigned int j = 0; j < width; j++)
		{
			re = 0;
			im = 0;
			for (unsigned int x = 0; x < height; x++)
			{
				for (unsigned int y = 0; y < width; y++)
				{
					temp = (float)(2 * PI * ((float)i * x / height + (float)j * y / width));
					re += (double)(img_gray[x][y] * cos(temp));
					im -= (double)(img_gray[x][y] * sin(temp));
				}
			}
			real[i][j] =(float)(re/(height*width));
			imaginary[i][j] = (float)(im/(height*width));
		}
	}
}

相应的IDFT公式可以转换为：
$$f(x,y)=\sum _{u=0}^{M-1}\sum _{v=0}^{N-1}F(u,v)[\cos 2\pi (ux/M+vy/N)+j\sin 2\pi (ux/M+vy/N)]$$
与上面的DFT算法类似，复杂度一样很高，对于尺寸较大的图像，即使在PC机上运行都需要很长的时间，对于一张250*130的图片，在我的电脑上单个函数的计算时间都达到了20秒。

//傅里叶逆变换，输入参数如DFT
void inverse_discrete_fourier_transform(unsigned char img_gray[][Width], float** real, float** imaginary, unsigned int height, unsigned int width)
{
	double re, temp;
	for (unsigned int i = 0; i < height; i++)
	{
		for (unsigned int j = 0; j < width; j++)
		{
			re = 0;
			for (unsigned int x = 0; x < height; x++)
			{
				for (unsigned int y = 0; y < width; y++)
				{
					temp = (float)(2 * PI * ((float)i*x/height+ (float)j*y/width));
					re += real[x][y] * cos( temp) - imaginary[x][y] * sin(temp);			
				}
			}
			img_gray[i][j] = (unsigned char)(re > 255 ? 255 : (re < 0 ? 0: re ));
		}
	}
}

算法的优化

上面虽然实现了DFT和IDFT，但是在单片机上完全不可以运行，因为用到了sin和cos函数，这两个函数在单片机上较难执行，所以需要一定的方法才能使用。在单片机上常见的sin和cos算法的实现通常是查表法，但是查表法会占用较多的空间，如果在精度上还有一定的要求并且数据比较随意，那么就会很难实现。
优化的数学原理是：
$$\begin{array}{r}\sin (a+b)=\sin (a)\cos (b)+\cos (a)\sin (b)\\ \cos (a)=\sin (a+90^{\circ })\\ \sin (a+180^{\circ })=-\sin (a)\end{array}$$
首先建立两个长度为10的数组，每隔10度存储对应的正弦值和每隔1度存储余弦值：

const float sin_table[] = {
    0.0,                                    //sin(0)
    0.17364817766693034885171662676931 ,    //sin(10)
    0.34202014332566873304409961468226 ,    //sin(20)
    0.5 ,                                   //sin(30)
    0.64278760968653932632264340990726 ,    //sin(40)
    0.76604444311897803520239265055542 ,    //sin(50)
    0.86602540378443864676372317075294 ,    //sin(60)
    0.93969262078590838405410927732473 ,    //sin(70)
    0.98480775301220805936674302458952 ,    //sin(80)
    1.0                                     //sin(90)
};

const float cos_table[] = {
    1.0 ,                                   //cos(0)
    0.99984769515639123915701155881391 ,    //cos(1)
    0.99939082701909573000624344004393 ,    //cos(2)
    0.99862953475457387378449205843944 ,    //cos(3)
    0.99756405025982424761316268064426 ,    //cos(4)
    0.99619469809174553229501040247389 ,    //cos(5)
    0.99452189536827333692269194498057 ,    //cos(6)
    0.99254615164132203498006158933058 ,    //cos(7)
    0.99026806874157031508377486734485 ,    //cos(8)
    0.98768834059513772619004024769344      //cos(9)
};

实现算法如下：

//弧度制下的角度转化成度数制下的角度时的比例因子,pi/180
const float hollyst = 0.017453292519943295769236907684886;
float qfsind(float x)
{
    int sig = 0;
    if (x > 0.0) {
        while (x >= 360.0) {
            x = x - 360.0*(unsigned int)(x/360);//因为计算的度数非常大，一次一次地减去360需要很久，
            //如果可以确定度数不大，可以直接如下所示
            //x = x - 360.0;
        }
    }
    else {
        while (x < 0.0) {
            x = x + 360.0 * (unsigned int)(x / 360);
        }
    }

    if (x >= 180.0) {
        sig = 1;
        x = x - 180.0;
    }

    x = (x > 90.0) ? (180.0 - x) : x;

    int a = x * 0.1;
    float b = x - 10 * a;
    float y = sin_table[a] * cos_table[(int)b] + b * hollyst * sin_table[9 - a];

    return (sig > 0) ? -y : y;
}
float qfcosd(float x)
{
    return qfsind(x + 90.0);
}

需要注意的是，上面使用的math.h函数使用的弧度制，即传入的是弧度，并且公式使用的也是弧度，而这两个优化函数实现的是角度制的计算，因此需要加上转换才能使用。
这两个函数已经经过测试，准确度还是不错的，与原本的函数对比，最大误差都是小数点后三位的，大多数情况下的误差都是小数点后五六位的，除非对精度要求非常非常高，否则完全可以替换使用。

//求取dft的实部和虚部
void discrete_fourier_transform(unsigned char img_gray[][Width], float **real, float **imaginary, unsigned int height, unsigned int width)
{
	double re, im, temp;
	for (unsigned int i = 0; i < height; i++)
	{
		for (unsigned int j = 0; j < width; j++)
		{
			re = 0;
			im = 0;
			for (unsigned int x = 0; x < height; x++)
			{
				for (unsigned int y = 0; y < width; y++)
				{
                    //不同之处
					temp = (float)(360 * ((float)i * x / height + (float)j * y / width));
                    re += (double)(img_gray[x][y] * qfcosd(temp));
					im -= (double)(img_gray[x][y] * qfsind(temp));
				}
			}
			real[i][j] =(float)( re/(height*width));
			imaginary[i][j] = (float)(im/(height*width));
		}
	}
}
//傅里叶逆变换，输入参数如DFT
void inverse_discrete_fourier_transform(unsigned char img_gray[][Width], float** real, float** imaginary, unsigned int height, unsigned int width)
{
	double re, temp;
	for (unsigned int i = 0; i < height; i++)
	{
		for (unsigned int j = 0; j < width; j++)
		{
			re = 0;
			for (unsigned int x = 0; x < height; x++)
			{
				for (unsigned int y = 0; y < width; y++)
				{
                    //不同之处
					temp = (float)(360 * ((float)i*x/height+ (float)j*y/width));
					re += real[x][y] * qfcosd( temp) - imaginary[x][y] * qfsind(temp);			
				}
			}
			img_gray[i][j] = (unsigned char)(re > 255 ? 255 : (re < 0 ? 0: re ));
		}
	}
}

编写一个测试函数，这个测试函数借用了malloc.h，一般有外部SDRAM的单片机也可以实现自己的malloc函数和free函数，故不作优化。

void test_dft(unsigned char img_gray[][Width], unsigned char result[][Width],unsigned int height,unsigned int width)
{
	float**real=(float**)malloc(sizeof(float*)*height);
	float**imaginary = (float**)malloc(sizeof(float*) * height);
	for (unsigned int i = 0; i < height; i++)
	{
		real[i] = (float*)malloc(sizeof(float) * width);
		imaginary[i] = (float*)malloc(sizeof(float) * width);
	}
	discrete_fourier_transform(img_gray, real, imaginary, Height, Width);
	inverse_discrete_fourier_transform(result , real, imaginary, Height, Width);
	for (unsigned int i = 0; i < height; i++)
	{
		free(real[i]);
		free(imaginary[i]);
	}
}

主函数部分延续上几篇文章的框架：

int main()
{
	Mat image_gray = imread("./test1.jpg", 0);
	unsigned char img_gray[Height][Width];
	unsigned char result[Height][Width];
	Mat2array(image_gray, img_gray, Height, Width);
	test_dft(img_gray, result,Height,Width);
	Mat img_test = array2Mat(result, Height, Width);
	Mat img_test2 = array2Mat(img_gray, Height, Width);
	imshow("image_test", img_test);
	imshow("image_test1", img_test2);
	waitKey(0);
	return 0;
}

测试结果

其中一张图片是输入前的灰度图，一张是经过DFT和IDFT的还原图，两者基本一致，肉眼看不出区别。
为什么不求取幅度谱来查看呢？
因为幅度谱的显示方法比较特殊，且为了突出显示效果，都会经过增减一定的数值、取对数、归一化等一系列操作（OpenCV官方例程），个人觉得没有太大的比较价值了。

附录

额外加上一个幅度谱的求取。首先搞清楚幅度谱的公式
$$|F(u,v)|=\sqrt{R(u,v{)}^2+I(u,v{)}^2}$$
即实部的平方加上虚部的平方后开方，对于单片机来说，平方不难实现，但是开方的难度较大，实现思路有两个：一是对于ARM内核的单片机，如STM32H7系列，其内部有专门的DSP算法实现，速度很快。二是使用卡马克算法进行开方，其精度还行，也是小数点后三四位，对于幅度在0-255的灰度而言，足够了。

//卡马克算法进行快速开平方
float Fsqrt(float x)
{
    int i;
    float x2, y;
    const float threehalfs = 1.5F;
    x2 = x * 0.5F;
    y = x;
    i = *(int*)&y;
    i = 0x5f375a86 - (i >> 1);
    y = *(float*)&i;
    y = y * (threehalfs - (x2 * y * y));
    y = y * (threehalfs - (x2 * y * y));
    y = y * (threehalfs - (x2 * y * y));
    return x * y;
}

因此在单片机中可以借助这个函数求取幅度谱。
快速开方算法

总结

这次使用的DFT和IDFT公式参考了贾永红的《数字图像处理》，因为在网上看到了一些与之不同的公式，故说明一下。
这仅仅是对DFT和IDFT基本公式的C语言转换而已，没有经过FFT的优化，所以更多地只是一次练手，并没有实际应用价值，但是本次查资料的过程中看到了挺不错的快速开方算法和快速sin函数实现，这两种优化方法更值得学习。