在數字信號和數字圖像領域，對頻域的研究是一個重要分支。

我們日常“加工”的圖像都是像素級，被稱為是圖像的空域數據。空域數據表征我們“可讀”的細節。如果我們將同一張圖像視為信號，進行頻譜分析，可以得到圖像的頻域數據。觀察下面這組圖，頻域圖中的亮點為低頻信號，代表圖像的大部分能量，也就是圖像的主體信息。暗點為高頻信號，代表圖像的邊緣和噪聲。從組圖可以看出，Degraded Goofy 與 Goofy 相比，近似的低頻信號保留住了 Goofy 的“輪廓”，而其高頻信號的增加使得背景噪點更加明顯。頻域分析使我們可以了解圖像的組成，進而做更多的抽象分析和細節處理。

Goofy and Degraded Goofy

實現圖像空域和頻域轉換的工具，就是傅立葉變換。由于圖像數據在空間上是離散的，我們使用傅立葉變換的離散形式 DFT（Discrete Fourier Transform）及其逆變換 IDFT（Inverse Discrete Fourier Transform)。Cooley-Tuckey在DFT的基礎上，開發了更快的算法 FFT（Fast Fourier Transform）。

DFT/FFT在數字圖像領域還有一些延伸應用。比如基于 DFT 的 DCT（Discrete Cosine Transform，離散余弦變換）就用在了圖像壓縮JPEG算法和圖像水印算法。

JPEG 編碼是通過色彩空間轉換、抽樣分塊、DCT 變換、量化編碼實現的。其中 DCT 變換的使用將圖像低頻信息和高頻信息區分開，在量化編碼過程中壓縮了少量低頻信息、大量高頻信息從而獲得尺寸上壓縮。從貓臉圖上可看出隨著壓縮比增大畫質會變差，但是主體信息還是得以保留。

圖像水印算法通過 DCT 將原圖轉換至頻域，選取合適的位置嵌入水印圖像信息，并通過 IDCT 轉換回原圖。這樣對原圖像的改變較小不易察覺，且水印通過操作可以被提取。

DFT/FFT 在深度學習領域也有延伸應用。比如利用 FFT 可以降低卷積計算量的特點，FFT_Conv 算法也成為常見的深度學習卷積算法。本文我們就來探究一下頻域算法的原理和優化策略。

DFT的原理及優化

公式

無論是多維的 DFT 運算，還是有基于 DFT 的 DCT/FFT_Conv, 底層的計算單元都是 DFT_1D。因此，DFT_1D 的優化是整個FFT類算子優化的基礎。

DFT_1D 的計算公式： 其中 為長度為 N 的輸入信號， 是 1 的 N 次根， 為長度為N 的輸出信號。

該公式的矩陣形式為：

單位復根的性質

DFT_1D 中的 是 1 的單位復根。直觀地看，就是將復平面劃分為 N 份，根據 k * n 的值逆時針掃過復平面的圓周。

單位復根有著周期性和對稱性，我們依據這兩個性質可以對W矩陣做大量的簡化，構成 DFT_1D 的快速算法的基礎。

周期性：

對稱性：

Cooley-Tuckey FFT算法

DFT_1D 的多種快速算法中，使用最頻繁的是 Cooley-Tuckey FFT 算法。算法采用用分治的思想，將輸入尺寸為 N 的序列，按照不同的基 radix，分解為 N/radix 個子序列，并對每個子序列再劃分，直到不能再被劃分為止。每一次劃分都可以得到一級 stage，將所有的級自下而上組合在一起，計算得到最后的輸出序列。

這里以 N = 8， radix=2 為例展示推理過程。

其中 為N=8 的序列， 為 DFT 輸出序列.

根據 DFT 的計算公式：

根據奇偶項拆開，分成兩個長度為 4 的序列 。

的 DFT 結果 乘以對應的旋轉因子 ，進行簡單的加減運算可以得到輸出 。

同理, 對 也做一樣的迭代， 都是 N=2 的序列，用他們的 DFT 結果進行組合運算可以得到 。

計算 N=2 的序列 , 因為 ，旋轉因子 。只要進行加減運算得到結果。

用算法圖形表示，每一層的計算會產生多個蝶形，因此該算法又被稱為蝶形算法。

這里我們要介紹碟形網絡的基本組成，對下文的分析有所幫助。

N=8 碟形算法圖

N=8 的計算序列被分成了 3 級，每一級 (stage) 有一個或多個塊 (section)，每個塊中包含了一個或者多個蝶形（butterfly) , 蝶形的計算就是 DFT 運算的 kernel。

每一個 stage 的計算順序：

• 取輸入
• 乘以轉換因子
• for section_num, for butterfly_num，執行radixN_kernel
• 寫入輸出

看 N=8 的蝶形算法圖，stage = 1 時，運算被分成了 4 個 section，每個section的 butterfly_num = 1 。stage = 2 時， section_num = 2，butterfly_num = 2。stage = 3時， section_num = 1， butterfly_num = 4。

可以觀察到，從左到右過程中 section_num 不斷減少， butterfly_num 不斷增加，蝶形群在“變大變密”，然而每一級總的碟形次數是不變的。

實際上，對于長度為 N ，radix = r 的算法，我們可以推得到：

為當前的 ，* sec/butterfly_stride 是每個section/butterfly* 的間隔。

這個算法可以將復雜度從 O(n^2) 下降到 O(nlogn)，顯得高效而優雅。我們基于蝶形算法，對于不同的radix進行算法的進一步劃分和優化，主要分為radix - 2 的冪次的和 radix – 非 2 的冪次兩類。

radix-2 的冪次優化

DFT_1D 的 kernel 即為矩陣形式中的 矩陣，我們對 radix_2^n的 kernel 進行分析。

背景里提到， DFT 公式的矩陣形式為：

其中 ~ 為乘以旋轉因子 后的輸入

當 radix = 2 時，由于 , radix_2 的 DFT 矩陣形式可以寫為：

當 radix = 4 時，由于 ，radix_4的DFT 矩陣形式可以寫為：

同理推得到 radix_8 的 kernel 為：

我們先來看訪存, 現代處理器對于計算性能的優化要優于對于訪存的優化, 在計算和訪存相近的場景下， 訪存通常是性能瓶頸。

DFT1D 中，對于不同基底的算法 r-2/r-4/r-8, 每一個 stage 有著相等的存取量: 2 * butterfly_num * radix = 2N, 而不同的基底對應的 stage 數有著明顯差異（ vs vs )。

因此對于 DFT , 在不顯著增加計算量的條件下， 選用較大的kernel會在訪存上取得明顯的優勢。觀察推導的 kernel 圖， r-2 的 kernel 每個蝶形對應 4次訪存操作和，2 次復數浮點加減運算。r-4 的 kernel 每個蝶形算法對應 8 次 load/store、8 次復數浮點加減操作（合并相同的運算），在計算量略增加的同時 stage 由 下降到 , 降低了總訪存的次數， 因此會有性能的提升。r-8 的 kerne l每個蝶形對應 16 次load/store、24 次復數浮點加法和8次浮點乘法。浮點乘法的存在使得計算代價有所上升， stage由 進一步下降到 ，但由于 N 日常并不會太大， r-4 到 r-8 的stage 減少不算明顯，所以優化有限

我們再來看計算的開銷. 減少計算的開銷通常有兩種辦法：減少多余的運算、并行化。

以 r-4 算法為例, kernel 部分的計算為：

• radix_4_first_stage(src, dst, sec_num, butterfly_num)
• radix_4_other_stage(src, dst, sec_num, butterfly_num)
• for Sec_num
• for butterfly_num
• raidx_4_kernel

radix4_first_stage 的數據由于 k=0, 旋轉因子都為 1 ，可以省去這部分復數乘法運算,單獨優化。radix4_other_stage 部分， 從第 2 個 stage 往后， butterfly_num = 4^(s-1) 都為 4 的倍數，而每個 butterfly 數組讀取/存儲都是間隔的。可以對最里層的循環做循環展開加向量化，實現 4 個或更多 butterfly 并行運算。循環展開和SIMD指令的使用不僅可以提高并行性， 也可以提升cacheline 利用的效率, 可以帶來較大的性能提升。以SM8150(armv8) 為例，r-4 的并行優化可以達到 r2 的 1.6x 的性能。

尺寸：1 * 2048（r2c) 環境：SM8150大核

總之，對于 radix-2^n 的優化，選用合適的 radix 以減少多 stage 帶來的訪存開銷，并且利用單位復根性質以及并行化降低計算的開銷，可以帶來較大的性能提升。

radix-非2的冪次優化

當輸入長度 N = radix1^m1 * radix2^m2… 且 radix 都不為 2 的冪次時，如果使用 naive的O(n^2) 算法， 性能就會急劇下降。常見的解決辦法對原長補 0、使用 radix_N 算法、特殊的 radix_N 算法(chirp-z transform)。補0至2的冪次方法對于大尺寸的輸入要增加很多運算量和存儲量， 而chirp-z transform 是用卷積計算DFT， 算法過于復雜。因此對非 2 的冪次radix-N 的優化也是必要的。

radix-N 計算流程和 radix-2 冪次一樣，我們同樣可以利用單位復根的周期性和對稱性，對 kernel 進行計算的簡化。以 radix-5 為例，radix-5 的DFT_kernel 為：

在復平面上根據x軸對稱，有相同的實部和相反的虛部。根據這個性質。如下圖所示，對于每一個 stage，可以合并公共項A，B，C，D，再根據公共項計算出該 stage 的輸出。

這種算法減少了很多重復的運算。同時，在 stage>=2 的時候，同樣對 butterfly 做循環展開加并行化，進一步減少計算的開銷。

radix-5 的優化思想可以外推至 radix-N 。對于 radix_N 的每一個 stage,計算流程為：

• 取輸入
• 乘以對應的轉換因子
• 計算公共項， radix_N 有 N-1個公共項
• 執行并行化的 radix_N_kernel
• 寫入輸出

其他優化

上述兩個章節描述的是 DFT_1D 的通用優化，在此基礎上還可以做更細致的優化，可以參考本文引用的論文。

• 對于全實數輸入的，由于輸入的虛部為 0， 進行旋轉因子以及radix_N_kernel 的復數運算時會有多余的運算和多余的存儲， 可以利用 split r2c 算法， 視為長度為 N/2 的復數序列， 計算 DFT 結果并進行 split操作得到 N 長實數序列的結果。
• 對于 radix-2 的冪次算法， 重新計算每個 stage 的輸入/輸出 stride 以取消第一級的位元翻轉可以進一步減少訪存的開銷。
• 對于 radix-N 算法， 在混合基框架下 N = radix1^m1 * radix2^m2， 合并較小的 radix 為大的 radix 以減少 stage。

DFT 延展算法的原理及優化

DCT 和FFT_conv 兩個典型的基于 DFT 延展的算法，DFT_1D/2D 的優化可以很好的用在這類算法中。

DCT

DCT算法（Discrete Cosine Transform， 離散余弦變換）可以看作是 DFT 取其正弦分量并經過工業校正的算法。DFT_1D 的計算公式為：

該算法naive實現是 O(n^2) 的，而我們將其轉換成 DFT_1D 算法，可以將算法復雜度降至 O(nlogn )。

基于 DFT 的 DCT 算法流程為：

• 對于 DCT 的輸入序列 x[n], 創建長為 2N 的輸入序列 y[n] 滿足 y[n] = x[n] + x[2N-n-1], 即做一個鏡像對稱。
• 對輸入序列 y[n] 進行 DFT 運算，得到輸出序列 Y[K]。
• 由 Y[K] 計算得到原輸入序列的輸出 X[K] 。

我們嘗試推導一下這個算法：

對 y[n] 依照 DFT 公式展開，整理展開的兩項并提取公共項 , 根據歐拉公式和誘導函數，整理非公共項 。可以看出得到的結果正是 x[k] 和與 k 有關的系數的乘積。這樣就可以通過先計算 得到 x[n] 的 DCT 輸出 。

在理解算法的基礎上，我們對 DFT_1D 的優化可以完整地應用到 DCT 上。DCT_2D 的計算過程是依次對行、列做 DCT_1D, 我們用多線程對 DCT_1D 進行并行，可以進一步優化算法。

FFT_conv

Conv 是深度學習最常見的運算，計算conv常用的方法有 IMG2COL+GEMM， Winograd， FFT_conv。三種算法都有各自的使用場景。

FFT_conv 的數學原理是時域中的循環卷積對應于其離散傅里葉變換的乘積. 如下圖所示， f 和 g 的卷積等同于將 f 和 g 各自做傅立葉變幻 F，進行點乘并通過傅立葉逆變換計算后的結果。

直觀的理論證明可下圖。

將卷積公式和離散傅立葉變換展開， 改變積分的順序并且替換變量， 可以證明結論。

注意這里的卷積是循環卷積， 和我們深度學習中常用的線性卷積是有區別的。利用循環卷積計算線性卷積的條件為循環卷積長度 L?| f |+| g |?1。因此我們要對 Feature Map 和 Kernel做zero-padding，并從最終結果中取有效的線性計算結果。

FFT_conv 算法的流程：

• 將 Feature Map 和 Kernel 都 zero-pad 到同一個尺寸,進行 DFT 轉換。
• 矩陣點乘
• 將計算結果通過 IDFT 計算出結果。

該算法將卷積轉換成點乘， 算法復雜度是 O(nlogn),小于卷積的 O(n^2), 在輸入的尺寸比較大時可以減少運算量，適用于大 kernel 的 conv 算法。

深度學習計算中， Kernel 的尺寸要遠小于 Feature Map， 因此 FFT_conv第一步的 zero-padding 會有很大的開銷，參考論文2里提到可以通過對 Feature map進行分塊， 分塊后的 Feature Map 和 Kernel 需要 padding 到的尺寸較小，可以大幅減小這一部分的開銷。優化后 fft_conv 的計算流程為：

• 合理安排緩存計算出合適的tile尺寸，對原圖進行分塊
• 分塊后的小圖和 kernel 進行 zero-padding , 并進行 DFT 運算
• 小圖矩陣點乘
• 進行逆運算并組合成大圖。

同時我們可以觀察到，FFT_conv 的核心計算模塊還是針對小圖的 DFT 運算， 因此我們可以將前一章節對 DFT 的優化代入此處，輔以多線程，進一步提升 FFT_Conv 的計算效率。