TMS320C6455二维FFT和IFFT实现

参考资料：

Rafael C. Gonzalez, Richard E. Woods, 数字图像处理（第三版）4.11，配套书内资源
Alan V.Oppenheim, Ronald W.Schafer, 离散时间信号处理（第三版）9.5
SPRUEB8B - TMS320C64x+ DSP Little-Endian DSP Library Programmer’s Reference
CCS Doc, 7.8 Image Analyzer

FFT原理简介

FFT的计算在任意一本《数字信号处理》的书中肯定都有讲，我看的是《离散时间信号处理》这本书。整个第9章对DFT和FFT做了详细的介绍。
最开始讲的Goertzel算法，引入一个数字序列，通过这一序列的递推公式求第N个值，等价于计算相应的DFT的结果。虽然这一算法的计算量仍然正比于 $N^2$ ,但也在一定程度上降低了计算复杂度。
而后的按时间抽取和按频率抽取的FFT算法在后面又归结为N为复合数的更为一般的FFT算法。按时间抽取的算法是比较好理解的。一个数字序列，根据标号的奇偶性分为两个序列，分别计算DFT，再由它们的结果计算整个序列的DFT就非常方便。其中可以灵活运用旋转因子的对称性和周期性，最后表示成一个“蝶形运算”。

FFT的计算中，数据存储也是一个需要审慎考虑的问题。蝶形运算两端的数据索引如果不需要发生改变，那样就可以把计算结果再存回原来的位置，而不需要额外开辟空间，这就被称为“同址计算（in-place）”。经过观察可以发现，如果输入是到位序的，然后用同址计算得到的结果就是正序的；如果输入是正序的，再用同址计算得到的结果就是到位序的。也有那种输入输出都是正序的算法，那样就没有采用同址计算。
我们熟悉的FFT算法，把序列分成奇偶序列递归地运算，就要求序列长度是2的幂次，这种算法称为“radix-2”算法；而如果按照类似的方法，把整个序列分成4组，即x[4n]、x[4n+1]、x[4n+2]、x[4n+3]，这被称为“radix-4”算法。那样所需要的级数（ $log_4N$ ）就更少，同时也要求序列的的长度是4的幂次。当然，还有将两种算法结合，主要采用radix-4，而在最后根据情况采用radix-2或者radix-4。

DSPLib配置

DSP库安装在哪也都无所谓，在Window->Preferences->Code Composer Studio->Products页面下可以找到安装的DSP库。
但是这个Product我不知道怎么用，就是我在实际的工程里即使勾选这个product好像也没啥用。
在跑例程的时候，我看到源文件中只包含了dsplib.h，所以只要包含它的路径和对应的静态库就行。
可以在Window->Preferences->Code Composer Studio->Build->Environment页面下添加DSP的安装路径。方便在工程中设置相应的包含路径和包含的静态库。

图像数据生成

用opencv-python读图像数据，调整图像分辨率，长宽同比例缩放，并在空白区域填充“0”，用于满足FFT对2的整数次幂的要求。生成输出.h文件和调整分辨率后的图像。

import cv2

img = cv2.imread("Fig0431(d)(blown_ic_crop).tif",cv2.IMREAD_GRAYSCALE)
fp_code = open('..\\FFT2\\image.h','w') #the path should adjust as your need

## the destinate image size
dst_height = 1024
dst_width = 1024

## the image size of raw image
height = img.shape[0]
width = img.shape[1]

#calculate the scale
scale = min(dst_height/height,dst_width/width)

#resize the image, using the same scale for both width and height
res = cv2.resize(img, None, fx=scale, fy=scale,interpolation=cv2.INTER_AREA)

#fill the reigon where is blank with black
res = cv2.copyMakeBorder(res,0,round(dst_height-height*scale), 0, round(dst_width-width*scale),cv2.BORDER_CONSTANT,value=[0])

## write data to .h file
linestride = int(dst_width/4)
fp_code.write('#pragma DATA_ALIGN(img_src, 8);\n')
fp_code.write('int img_src[] = {\n    ')
for i in range(dst_height):
    for k in range(4):
        for j in range(linestride):
            if(i==dst_height-1 and k == 3 and j == linestride-1):
                fp_code.write(str(res[i][k*linestride+j]) + ', 0')
            else:
                fp_code.write(str(res[i][k*linestride+j]) + ', 0, ')
        fp_code.write('\n    ')
fp_code.write('};\n')

## write image
cv2.imwrite("img.png",res)

DSPLib中的FFT和IFFT

DSPLib中的函数可以从源文件中的注释进行理解，也可以直接看SPRUEB8B这个文档。FFT根据输入和输出数据的位宽也分好多种，我实现的是那种输入输出都是32位的fft，没有缩放系数，也就是在"dsplib安装路径\packages\ti\dsplib\src"路径下的DSP_fft32x32和DSP_ifft32x32。
在上面的路径下有四个测试例程的文件夹还有三种FFT的实现方式，分别是纯C语言实现（_cn）、带有编译器指令的C语言实现（_i）和汇编实现（_sa）。这三种方式中，直接采用汇编实现的代码效率最高，速度最快。而DSP库默认调用的也就是这个汇编的FFT。
这三种实现方式，采用的计算方式都是相同的，只是采用专门的汇编指令有更高的效率。所以想要理解这个代码可以直接看_cn版本的。蝶形运算的示意图大概是下面这个样子：

输入输出都是正序，没有同址计算。采用radix-4的算法，一个蝶形运算需要4个数。在第一级运算中，这四个数的标号间隔为N/4（类比图里radix-2是N/2），这在代码中用变量stride表示，后续每级逐级减小。
IFFT的计算和FFT非常像，我们可以用类似的方式计算IFFT，《数字图像处理》的4.11.2给出了计算方法。

$N{x^*}[n] = \sum\limits_{n = 0}^{N - 1} { {X^*}[k]{e^{ - j2\pi kn/N} } }$

一般情况下 $x^*[n]=x[n]$ ，所以我们只需要对原复数序列取共轭后做DFT，再把得到的结果除以“N”就能得到IFFT的结果。

二维FFT和IFFT实现

在《数字图像处理》4.11.1中介绍了二维DFT的可分性，即要实现二维DFT，可以对行列两个维度分别做DFT。在代码中实现时，可以先对行做FFT，再将结果转置，之后再对行做FFT，最后再转置回去就能够得到对应的二维DFT的结果。\

/*
 * @func fft2, because of transpose function, this function can be only used when width is equal to height!!!
 *
 * @param w = twiddle factors
 * @param psrc = source pointer, point to the image data to be processed
 * @param ptmp = a pointer point to a data space, witch size is as same as the source image, just for temporary use
 * @param width = image width
 * @param height = image height
 */
void DSP_fft2(int *w, int *psrc, int *ptmp, int width, int height){
    int *ps = psrc;
    int *pd = ptmp;
    int i;

    // fft for the row
    for(i=0;i<height;i++){
        DSP_fft32x32(w, width, ps, pd);
        ps = ps + 2*width;
        pd = pd + 2*width;
    }

    //transpose
    transpose(ptmp, width, height);

    //fft for the column
    ps = ptmp;
    pd = psrc;
    for(i=0;i<width;i++){
        DSP_fft32x32(w, width, ps, pd);
        ps = ps + 2*width;
        pd = pd + 2*width;
    }

    //transpose
    transpose(psrc, width, height);
}

这里的转置就是一般的转置，不是共轭转置。如果图像的行和列长度是相等的，转置只需要交换行列坐标相对的两个位置的数据。而如果图像的行和列的长度不相等，则需要额外开辟空间用来存放转置的数据，额外的存储开销较大，因此我就没写这部分代码，而是要求图像行列等长。
DSP库中提供的IFFT仅是在FFT的基础上进行了一点修改，改了一些符号，实现了取复共轭的计算，但是没有除以序列长度。因此在每次调用DSP_ifft32x32后需要对结果再除以“N”才是IFFT的结果。

/*
 * @func ifft2, because of transpose function, this function can be only used when width is equal to height!!!
 *
 * @param w  = twiddle factors
 * @param psrc = source pointer, point to the image data to be processed
 * @param ptmp = a pointer point to a data space, witch size is as same as the source image, just for temporary use
 * @param width = image width
 * @param height = image height
 */
void DSP_ifft2(int *w, int *psrc, int *ptmp, int width, int height){
    int *ps = psrc;
    int *pd = ptmp;
    int i;

    //ifft for the row
    for(i=0;i<height;i++){
        DSP_ifft32x32(w, width, ps, pd);
        ps = ps + 2*width;
        pd = pd + 2*width;
    }
    for(i=0;i<2*width*height;i++){
        ptmp[i] = ptmp[i]/width;
    }

    //transpose
    transpose(ptmp, width, height);

    //ifft for the column
    ps = ptmp;
    pd = psrc;
    for(i=0;i<width;i++){
        DSP_ifft32x32(w, width, ps, pd);
        ps = ps + 2*width;
        pd = pd + 2*width;
    }
    for(i=0;i<2*width*height;i++){
        psrc[i] = psrc[i]/width;
    }

    //transpose
    transpose(psrc, width, height);
}