目录
参考资料：
- 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
N^2
N2,但也在一定程度上降低了计算复杂度。
而后的按时间抽取和按频率抽取的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”算法。那样所需要的级数（
l
o
g
4
N
log_4N
log4N）就更少，同时也要求序列的的长度是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
]
=
∑
n
=
0
N
−
1
X
∗
[
k
]
e
−
j
2
π
k
n
/
N
N{x^*}[n] = \sum\limits_{n = 0}^{N - 1} {{X^*}[k]{e^{ - j2\pi kn/N}}}
Nx∗[n]=n=0∑N−1X∗[k]e−j2πkn/N
一般情况下
x
∗
[
n
]
=
x
[
n
]
x^*[n]=x[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);
}
图像分析工具（Image Analyzer）
CCS自带图像分析工具。在调试界面，从Tools->Image Analyzer打开，电梯Image窗口内的区域，可以编辑Properties窗口中的属性。
- Image Format: 输入图像格式，取决于具体的数据格式
- Number of pixels per line: 每行像素
- Number of lines: 图像行数
- Data format: 数据格式，packed就是一个像素的数据是连续存放的，还有一种planar的就是同一通道的数据是连续存放的
- Pixel stride (bytes): 像素步进，就是每个像素占几个字节
- Red/Green/Blue/Alpha mask: “1”有效，可以让ARGB从几个字节中找到对应位置的值，可以重叠，RGB三个通道重叠就是灰度图像。
- Line stride (bytes): 每行占多少字节
- Start address: 图像起始地址
- Read data as: 读数据的格式，根据图像的数据格式设置，int(32bit), short(16bit), char(8bit)
- Start address: 图像起始地址
- Read data as: 读数据的格式，根据图像的数据格式设置，int(32bit), short(16bit), char(8bit)
设置完成后在Image窗口中点刷新按钮，即可刷新显示图像。
测试结果
预先存入的图像：
经过FFT变换，实际数据已经超出灰度范围，无法显示正常的FFT频谱：
IFFT变换后的图像：