FFT原理簡介

FFT的計算在任意一本《數字信号處理》的書中肯定都有講，我看的是《離散時間信号處理》這本書。整個第9章對DFT和FFT做了詳細的介紹。

最開始講的Goertzel算法，引入一個數字序列，通過這一序列的遞推公式求第N個值，等價于計算相應的DFT的結果。雖然這一算法的計算量仍然正比于 N 2 N^2 N2,但也在一定程度上降低了計算複雜度。

而後的按時間抽取和按頻率抽取的FFT算法在後面又歸結為N為複合數的更為一般的FFT算法。按時間抽取的算法是比較好了解的。一個數字序列，根據标号的奇偶性分為兩個序列，分别計算DFT，再由它們的結果計算整個序列的DFT就非常友善。其中可以靈活運用旋轉因子的對稱性和周期性，最後表示成一個“蝶形運算”。

TMS320C6455二維FFT和IFFT實作FFT原理簡介DSPLib配置圖像資料生成DSPLib中的FFT和IFFT二維FFT和IFFT實作圖像分析工具（Image Analyzer）測試結果相關資源連結

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的安裝路徑。友善在工程中設定相應的包含路徑和包含的靜态庫。

TMS320C6455二維FFT和IFFT實作FFT原理簡介DSPLib配置圖像資料生成DSPLib中的FFT和IFFT二維FFT和IFFT實作圖像分析工具（Image Analyzer）測試結果相關資源連結

圖像資料生成

用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版本的。蝶形運算的示意圖大概是下面這個樣子：

TMS320C6455二維FFT和IFFT實作FFT原理簡介DSPLib配置圖像資料生成DSPLib中的FFT和IFFT二維FFT和IFFT實作圖像分析工具（Image Analyzer）測試結果相關資源連結

輸入輸出都是正序，沒有同址計算。采用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視窗中的屬性。

TMS320C6455二維FFT和IFFT實作FFT原理簡介DSPLib配置圖像資料生成DSPLib中的FFT和IFFT二維FFT和IFFT實作圖像分析工具（Image Analyzer）測試結果相關資源連結

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視窗中點重新整理按鈕，即可重新整理顯示圖像。

測試結果

預先存入的圖像：

TMS320C6455二維FFT和IFFT實作FFT原理簡介DSPLib配置圖像資料生成DSPLib中的FFT和IFFT二維FFT和IFFT實作圖像分析工具（Image Analyzer）測試結果相關資源連結

經過FFT變換，實際資料已經超出灰階範圍，無法顯示正常的FFT頻譜：

TMS320C6455二維FFT和IFFT實作FFT原理簡介DSPLib配置圖像資料生成DSPLib中的FFT和IFFT二維FFT和IFFT實作圖像分析工具（Image Analyzer）測試結果相關資源連結

IFFT變換後的圖像：

TMS320C6455二維FFT和IFFT實作FFT原理簡介DSPLib配置圖像資料生成DSPLib中的FFT和IFFT二維FFT和IFFT實作圖像分析工具（Image Analyzer）測試結果相關資源連結

TMS320C6455二維FFT和IFFT實作FFT原理簡介DSPLib配置圖像資料生成DSPLib中的FFT和IFFT二維FFT和IFFT實作圖像分析工具（Image Analyzer）測試結果相關資源連結

目錄

FFT原理簡介

DSPLib配置

圖像資料生成

DSPLib中的FFT和IFFT

二維FFT和IFFT實作

圖像分析工具（Image Analyzer）

測試結果

相關資源連結

繼續閱讀

檔案目标uboot的配置過程分析

UBoot的工作

藍橋杯單片機經曆

【藍橋杯單片機06】矩陣鍵盤的基本控制

藍橋杯單片機組學習總結（一）

Build大會沒有了Windows，微軟祭出哪些AI殺手锏？

“2014全球未來網絡與SDN技術大會”精彩剪輯（上）

轉：2008年微軟Windows硬體工程(WinHEC)大會

2012年全球最愚蠢的設計第一是微軟，第二還是微軟

2014年軟體業經濟運作情況

UART學習之路（一）基本概念

CAN(一)-基本概念CAN 簡介

彙編語言基本概念（續13）

The first led（led逐個點亮過程

無人機--飛控科普

STM32-UART序列槽通信STM32-UART序列槽通信