利用傅立叶变换编程做出下面图像的二维DFT（说明：该图像像素点个数：64 * 64，中央黑色方块的像素点个数：8 * 8。黑色部分的灰度值是0，白色部分的灰度值是255。）。
打印出该图像经过变换后的幅度谱图像。
一. 方法描述
1．    DFT算法步骤：
 1）设f(x,y)是原图像的函数，x和y分别表示第x行第y列。首先对f(x,y)进行中心化；
 2) 设F（u,v）是傅立叶变换后的函数，u和v分别表示第u行第v列，根据欧拉公式 进行变换，得
 F(u,v)=a+jb
2．FFT算法步骤：
   1）迭代总的次数： 
         其中，N代表序列的长度。
   2）求对偶结点： 
        其中，e代表迭代到第几次，k代表本结点在序列中的号数。
   3）求加权系数W，关键是求P。步骤如下：
        a.将k写成r位的二进制数。
        b.右移r-e位，左边补零。
        c.将右移的二进制数比特倒转，结果换算成十进制就是所要求的P值。
   4）重新排序
        a.将最后一次迭代的结果Xe(k)中的序号数k写成二进制数
        b.进行比特倒转
        c.转换为十进制，就得到与X(k)对应的X(m)的序号数
   5）主函数：
a.首先求出每一行的DFFT，并存储在中间矩阵数组内；
b.求出中间数组矩阵的每一列的FFT，得到的结果就是2DFFT。
c.求出每个像素点的幅值|H(u,v)|。
d.作对数变换：D(u,v)=lg(1+|H(u,v)|)
e.量化成0 ~ 255之间的值。
3．硬件环境及程序说明：
本程序使用Visual C++ 6.0编写，在测控实验室IBM兼容PC机Windows2000Server操作系统上面运行实现。为了使程序有好的可移植性，核心算法很少调用Windows提供的API函数，大部分使用ANSI C提供的标准库函数。
1）应用程序框架
    应用程序框架由Visual C++中提供的MFC自动生成。首先提供了两个类：应用程序类和CavgbmpDlg类。
2）DFT变换
    DFT变换使用CavgbmpDlg类中下面3个函数：
    void CAvgbmpDlg::OnDft()：响应开始做DFT变换按钮事件的函数。该函数首先声明了目标矩阵F，然后初始化原始矩阵的每个点。填入灰度值。其中白色部分填入255，黑色部分填入0。然后调用DFT(f)进行傅立叶变换。最后调用DrawDFT()函数将F的每个点的灰度值画到屏幕上。
    Den CAvgbmpDlg::DFT(BYTE s[64][64])：本函数是DFT变换的核心算法函数。参数s[][]是要变换的原始矩阵。首先对s进行中心化，然后根据公式进行变换，最后求模，得到变换后的结果。
    void CAvgbmpDlg::DrawDFT(int x, int y, Den aa)：本函数的功能是将变换后的图像画到屏幕上。参数x,y是屏幕的坐标，aa是目标矩阵（变换后的矩阵）。
3）FFT变换
FFT变换用到以下的函数：
void CAvgbmpDlg::OnFft()：响应开始做FFT变换按钮事件的函数。该函数首先对目标矩阵和原始矩阵进行初始化。然后调用kbFFT()函数对原始矩阵进行行FFT变换，然后对得到的中间矩阵作列FFT变换。最后对所得结果进行量化，将结果画到屏幕上。原理是：二维离散傅立叶变换可分离性。
void CAvgbmpDlg::kbFFT(double pr[N], double pi[N], int n, int k, double fr[N], double fi[N], int l, int il)：本函数是FFT变换的核心函数。参数说明如下：
pr[N]：双精度实型一维数组，长度为N，当l=0时，存放采样输入的实部，返回时，存放离散傅立叶变换的模；当l=1时，存放傅立叶变换的N个实部，返回时，存放逆傅立叶变换的模。
Pi[N]：双精度实型一维数组，长度为N，当l=0时，存放采样输入的虚部，返回时，存放离散傅立叶变换的幅角；当l=1时，存放傅立叶变换的N个虚部，返回时，存放逆傅立叶变换的幅角，单位为度。
n：整型变量，表示输入的点数。
k：整型变量，表示迭代次数，相当于公式中的r。
fr[N]：双精度实型一维数组，长度为N，当l=0时，返回傅立叶变换的实部；当l=1时，返回逆傅立叶变换的实部。
fi[N]：双精度实型一维数组，长度为N，当l=0时，返回傅立叶变换的虚部；当l=1时，返回逆傅立叶变换的虚部。
l：整型变量，若l=0，表示不要求本函数计算傅立叶变换；若l=1，表示不要求本函数计算逆傅立叶变换。
il：整型变量，若il=0，表示不要求本函数计算傅立叶变换或逆变换的模与幅角；若il=1，表示要求本函数计算傅立叶变换或逆变换的模与幅角。
4）剩余的问题
浮点运算以及数据类型强制转换可能带来一些问题，出现误差。
例如，开始时候计算迭代次数r，根据公式 ，写出程序：r=log(64) / log(2);计算结果是5，而正确结果是6，也就是说计算机需要迭代6次才能得到结果。经过断点调试跟踪才发现本行代码计算结果有误。于是不得不在KbFFT()函数中再加入参数r。
另外，DFT算出的F(u,v)值很小，需要放大10倍才可以出现比较正常的结果。FFT算出的F(u,v)很大，需要缩小10倍才可以出现比较正常的结果。原因还不详。
二．程序示例如下：
//源程序：
//--------------------------------------------------------
//核心数据结构
typedef struct{ //存储变换后的目标矩阵
 BYTE A[64][64];
}Den;
class CAvgbmpDlg : public CDialog  //对话框类。
{
// Construction
public:
 CDC* pDC; //获取标准输出设备
 CAvgbmpDlg(CWnd* pParent = NULL); // standard constructor
 Den DFT(BYTE s[64][64]);
 void DrawDFT(int x, int y,Den aa);
    void kbFFT(double pr[N], double pi[N], int n,int k,  double fr[N], double fi[N], int l, int il);
};
void CAvgbmpDlg::kbFFT(double&nbs