浅谈原理
在图像中像素点为一个个离散的点，因此可以利用离散信号x(nT)来计算。又因为对图像的二维快速傅立叶变化相当于对图像的像素矩阵进行横和竖两次一维的快速傅立叶变化，我们这里原理按照一维快速傅立叶变换说明。
首先，有限长离散信号x(n)，n=0，1，…，N-1的DFT定义为：
x(t)=∑n=0N−1x(n)WknN,(k=0,1,…,N−1,WN=e−j2∏n)x(t)=∑n=0N−1x(n)WNkn,(k=0,1,…,N−1,WN=e−j2∏n)
将x(n)分解为偶数与奇数的两个序列之和，即
x(n)=x1（n）+x2（n）x(n)=x1（n）+x2（n）
x1(n)和x2(n)的长度都是N／2，x1(n)是偶数序列，x2(n)是奇数序列，则
X(k)=∑n=0N2−1x1W2nkN+∑n=0N2x2W(2n+1)kN,(k=0,1,…,N−1)X(k)=∑n=0N2−1x1WN2nk+∑n=0N2x2WN(2n+1)k,(k=0,1,…,N−1)
由于
W2nkN=e−j2∏N2nk=e−j2∏N/2nk=WnkN/2WN2nk=e−j2∏N2nk=e−j2∏N/2nk=WN/2nk
则原式等于
X(k)=∑n=0N2−1x1(n)WnkN/2+WkN∑n=0N2−1x2(n)WnkN/2=X1(k)+WkNX2(k),(k=0,1,…,N−1)X(k)=∑n=0N2−1x1(n)WN/2nk+WNk∑n=0N2−1x2(n)WN/2nk=X1(k)+WNkX2(k),(k=0,1,…,N−1)
其中X1(k)和X2(k)分别为x1(n)和x2(n)的N／2点DFT。由于X1(k)和X2(k)均以N／2为周期，且WN k+N/2=-WN k，所以X(k)又可表示为：
X(k)=X1(k)+WkNX2(k),(k=0,1,…,N/2−1)X(k)=X1(k)+WNkX2(k),(k=0,1,…,N/2−1)
X(k+N2)=X1(k)−WkNX2(k),(k=0,1,…,N/2−1)X(k+N2)=X1(k)−WNkX2(k),(k=0,1,…,N/2−1)
上式的运算称之为蝶形运算，依此类推，经过m-1次分解，最后将N点DFT分解为N／2个两点DFT。下图为8点FFT的分解流程：
图像快速傅立叶变换用到的函数和类
从原理中可以看出，我们在对图像进行快速傅立叶变换时需要图像的长和宽都为2的幂次方，所以首先我们需要对图像的长和宽进行拓展，矩阵拓展的部分补0。
为方便的找到长和宽需要拓展的最小2幂值，将寻找最小2次幂值函数写成一个函数的形式，输入为int的值，输出也为int的值。
寻找最小2次幂值具体函数如下：
public int find2n(int num){double b = Math.Log(num, 2);int c = (int)(Math.Log(num, 2));if (b - c != 0){num = (int)Math.Pow(2, (int)c + 1);}return num;}
计算出需要拓展的长和宽后，我们需要对原图进行拓展补0。将拓展写成一个函数的形式，输入为byte[]格式的图像、图像的长和宽和需要拓展的长和宽，输出也为byte[]格式的图像。这里将图像写成一维的形式还是二维形式一样，对应关系为：[i, j]＝[i * imageWidth + j]。
拓展具体函数如下：
//imageWidth和imageHeight为原图的宽度和长度；Width和Height为需拓展的宽度和长度public byte[] expand(byte[] image, int imageWidth, int imageHeight, int Width, int Height){byte[,] matrix = new byte[Height, Width];for (int i = 0; i < imageHeight; i++){for (int j = 0; j < imageWidth; j++){matrix[i, j] = image[i * imageWidth + j];}}for (int i = imageHeight; i < Height; i++){for (int j = imageWidth; j < Width; j++){matrix[i, j] = 0;}}byte[] result = new byte[Width * Height];for (int i = 0; i < Height; i++){for (int j = 0; j < Width; j++){result[i * Width1 + j] = matrix[i, j];}}return result;}
又由于在计算中会涉及到频域计算，需要复数的计算，而C#中没有将复数设计为一种内部数据类型，所以先要自己设计并处理复数运算的C#类Complex。
在解决方案资源管理器中右键点击添加新建项，在“名称“内输入类的名称，我们命名为Complex.cs，然后单击“添加“完成Complex类的添加。
具体Complex类定义如下：
class Complex{//复数的实部private double real = 0.0;//复数的虚部private double imaginary = 0.0;//实部的属性public double Real{get{return real;}set{real = value;}}//虚部的属性public double Imaginary{get{return imaginary;}set{imaginary = value;}}//基本构造函数public Complex(){}//指定值得构造函数public Complex(double dbreal, double dbimag){real = dbreal;imaginary = dbimag;}//复制构造函数public Complex(Complex other){real = other.real;imaginary = other.imaginary;}//重载+运算符public static Complex operator +(Complex comp1, Complex comp2){return comp1.Add(comp2);}//重载-运算符public static Complex operator -(Complex compl, Complex comp2){return compl.Subtract(comp2);}//重载*运算符public static Complex operator *(Complex compl, Complex comp2){return compl.Multiply(comp2);}//实现复数加法public Complex Add(Complex comp){double x = real + comp.real;double y = imaginary + comp.imaginary;return new Complex(x, y);}//实现复数减法public Complex Subtract(Complex comp){double x = real - comp.real;double y = imaginary - comp.imaginary;return new Complex(x, y);}//实现复数乘法public Complex Multiply(Complex comp){double x = real * comp.real - imaginary * comp.imaginary;double y = real * comp.imaginary + imaginary * comp.real;return new Complex(x, y);}//求幅度public double Abs(){//取得实部的绝对值double x = Math.Abs(real);//取得虚部的绝对值double y = Math.Abs(imaginary);//实部为0if (real == 0){return y;}//虚部为0if (imaginary == 0){return x;}//计算模if (x > y){return (x * Math.Sqrt(1 + (y / x) * (y / x)));}else{return (y * Math.Sqrt(1 + (y / x) * (y / x)));}}//求相位角public double Angle(){//实数和虚数都为0if (real == 0 && imaginary == 0)return 0;if (real == 0){//实部位0if (imaginary > 0)return Math.PI / 2;elsereturn -Math.PI / 2;}else{if (real > 0){//实部大于0return Math.Atan2(imaginary, real);}else{//实部小于0if (imaginary >= 0)return Math.Atan2(imaginary, real) + Math.PI;elsereturn Math.Atan2(imaginary, real) - Math.PI;}}}//共轭复数public Complex conjugate(){return new Complex(this.real, -this.imaginary);}}
C#实现图像快速傅立叶变换
我们将图像快速傅立叶变换也写成一个函数的形式，输入为byte[]格式的图像和图像的长和宽，输出为Complex[]格式的频域信息。
图像快速傅立叶变换具体函数如下：
//图像傅里叶变换private Complex[] FFT(byte[] imageData, int imageWidth, int imageHeight){int bytes = imageWidth * imageHeight;byte[] bmpValues = new byte[bytes];Complex[] tempCom1 = new Complex[bytes];bmpValues = (byte[])imageData.Clone();//赋值：把实数变为复数，即虚部为0for (int i = 0; i < bytes; i++){tempCom1[i] = new Complex(bmpValues[i], 0);}//水平方向快速傅里叶变换Complex[] tempCom2 = new Complex[imageWidth];Complex[] tempCom3 = new Complex[imageWidth];for (int i = 0; i < imageHeight; i++){for (int j = 0; j < imageWidth; j++){tempCom2[j] = tempCom1[i * imageWidth + j];}//调用一维傅里叶变换tempCom3 = fft(tempCom2, imageWidth);Complex[] tempCom9 = new Complex[imageWidth];for (int j = 0; j < imageWidth; j++){tempCom9[j] = tempCom3[j];}//将结果赋值回去for (int j = 0; j < imageWidth; j++){tempCom1[i * imageWidth + j] = tempCom3[j];}}//垂直方向傅里叶变换Complex[] tempCom4 = new Complex[imageHeight];Complex[] tempCom5 = new Complex[imageHeight];for (int i = 0; i < imageWidth; i++){for (int j = 0; j < imageHeight; j++){tempCom4[j] = tempCom1[j * imageWidth + i];}//调用一维傅里叶变换tempCom5 = fft(tempCom4, imageHeight);//把结果赋值回去for (int j = 0; j < imageHeight; j++){tempCom1[j * imageWidth + i] = tempCom5[j];}}return tempCom1;}
在图像傅立叶变换时中调用的一维快速傅立叶变换的函数，其输入为Complex[]格式的频域信息，输出也为Complex[]格式的频域信息。
一维快速傅立叶变换具体函数如下：
//一维傅里叶变换private Complex[] fft(Complex[] sourceData, int countN){//fft的级数int r = Convert.ToInt32(Math.Log(countN, 2));Complex[] w = new Complex[countN / 2];Complex[] interVar1 = new Complex[countN];Complex[] interVar2 = new Complex[countN];interVar1 = (Complex[])sourceData.Clone();//求加权系数wfor (int i = 0; i < countN / 2; i++){double angle = -i * Math.PI * 2 / countN;w[i] = new Complex(Math.Cos(angle), Math.Sin(angle));}//蝶形运算for (int i = 0; i < r; i++){int interval = 1 << i;int halfN = 1 << (r - i);//对每级的每一组点循环for (int j = 0; j < interval; j++){int gap = j * halfN;//对每级的每一点循环for (uint k = 0; k < halfN / 2; k++){//运行蝶形算法interVar2[k + gap] = interVar1[k + gap] + interVar1[k + gap + halfN / 2];interVar2[k + gap + halfN / 2] = (interVar1[k + gap] - interVar1[k + gap + halfN / 2]) * w[k * interval];}}interVar1 = (Complex[])interVar2.Clone();}//按位取反for (uint j = 0; j < countN; j++){uint rev = 0;uint num = j;//重新排序for (uint i = 0; i < r; i++){rev <<= 1;rev |= num & 1;num >>= 1;}interVar2[rev] = interVar1[j];}return interVar2;}
C#实现图像快速傅立叶反变换
相同的，我们将图像快速傅立叶反变换也写成一个函数的形式，输入为Complex[]格式的频域信息和图像的长和宽，输出为double[]格式的图像。这里输出double[]而不是byte[]是因为有时候后续图像处理操作有较高的精度要求。
图像快速傅立叶反变换具体函数如下：
//傅里叶反变换private double[] IFFT(Complex[] freData, int imageWidth, int imageHeight){int bytes = imageWidth * imageHeight;byte[] bmpValues = new byte[bytes];Complex[] tempCom1 = new Complex[bytes];tempCom1 = (Complex[])freData.Clone();//水平方向傅里叶逆变化Complex[] tempCom2 = new Complex[imageWidth];Complex[] tempCom3 = new Complex[imageWidth];for (int i = 0; i < imageHeight; i++){//得到水平方向复数序列{for (int j = 0; j < imageWidth; j++){tempCom2[j] = tempCom1[i * imageWidth + j];}//调用一维傅里叶变换tempCom3 = ifft(tempCom2, imageWidth);//把结果赋值回去for (int j = 0; j < imageWidth; j++){tempCom1[i * imageWidth + j] = tempCom3[j];}}}//垂直方向傅里叶逆变换Complex[] tempCom4 = new Complex[imageHeight];Complex[] tempCom5 = new Complex[imageHeight];for (int i = 0; i < imageWidth; i++){//得到垂直方向复数序列for (int j = 0; j < imageHeight; j++){tempCom4[j] = tempCom1[j * imageWidth + i];}//调用一维傅里叶变换tempCom5 = ifft(tempCom4, imageHeight);//把结果赋值回去for (int j = 0; j < imageHeight; j++){tempCom1[j * imageWidth + i] = tempCom5[j];}}//赋值：把复数转换为实数，只保留复数的实数部分double[] tempDouble = new double[bytes];for (int i = 0; i < bytes; i++){tempDouble[i] = tempCom1[i].Real;}return tempDouble;}
在图像傅立叶反变换时中调用的一维快速傅立叶反变换的函数，其输入为Complex[]格式的频域信息，输出也为Complex[]格式的频域信息。
一维快速傅立叶反变换具体函数如下：
//一维傅里叶反变换private Complex[] ifft(Complex[] sourceData, int countN){//共轭变换for (uint i = 0; i < countN; i++){sourceData[i] = sourceData[i].conjugate();}Complex[] interVar = new Complex[countN];//调用快速傅里叶变换interVar = fft(sourceData, countN);//共轭变换，并除以长度for (uint i = 0; i < countN; i++){interVar[i] = new Complex(interVar[i].Real / countN, -interVar[i].Imaginary / countN);}return interVar;}
