在數(shù)字信號和數(shù)字圖像領(lǐng)域，對頻域的研究是一個(gè)重要分支。

我們?nèi)粘！凹庸ぁ钡膱D像都是像素級，被稱為是圖像的空域數(shù)據(jù)?？沼驍?shù)據(jù)表征我們“可讀”的細(xì)節(jié)。如果我們將同一張圖像視為信號，進(jìn)行頻譜分析，可以得到圖像的頻域數(shù)據(jù)。觀察下面這組圖，頻域圖中的亮點(diǎn)為低頻信號，代表圖像的大部分能量，也就是圖像的主體信息。暗點(diǎn)為高頻信號，代表圖像的邊緣和噪聲。從組圖可以看出，Degraded Goofy 與 Goofy 相比，近似的低頻信號保留住了 Goofy 的“輪廓”，而其高頻信號的增加使得背景噪點(diǎn)更加明顯。頻域分析使我們可以了解圖像的組成，進(jìn)而做更多的抽象分析和細(xì)節(jié)處理。

Goofy and Degraded Goofy

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

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

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

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

DFT/FFT 在深度學(xué)習(xí)領(lǐng)域也有延伸應(yīng)用。比如利用 FFT 可以降低卷積計(jì)算量的特點(diǎn)，F(xiàn)FT_Conv 算法也成為常見的深度學(xué)習(xí)卷積算法。本文我們就來探究一下頻域算法的原理和優(yōu)化策略。

DFT的原理及優(yōu)化

公式

無論是多維的 DFT 運(yùn)算，還是有基于 DFT 的 DCT/FFT_Conv, 底層的計(jì)算單元都是 DFT_1D。因此，DFT_1D 的優(yōu)化是整個(gè)FFT類算子優(yōu)化的基礎(chǔ)。

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

該公式的矩陣形式為：

單位復(fù)根的性質(zhì)

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

單位復(fù)根有著周期性和對稱性，我們依據(jù)這兩個(gè)性質(zhì)可以對W矩陣做大量的簡化，構(gòu)成 DFT_1D 的快速算法的基礎(chǔ)。

周期性：

對稱性：

Cooley-Tuckey FFT算法

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

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

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

根據(jù) DFT 的計(jì)算公式：

根據(jù)奇偶項(xiàng)拆開，分成兩個(gè)長度為 4 的序列 。

的 DFT 結(jié)果 乘以對應(yīng)的旋轉(zhuǎn)因子 ，進(jìn)行簡單的加減運(yùn)算可以得到輸出 。

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

計(jì)算 N=2 的序列 , 因?yàn)? ，旋轉(zhuǎn)因子 。只要進(jìn)行加減運(yùn)算得到結(jié)果。

用算法圖形表示，每一層的計(jì)算會產(chǎn)生多個(gè)蝶形，因此該算法又被稱為蝶形算法。

這里我們要介紹碟形網(wǎng)絡(luò)的基本組成，對下文的分析有所幫助。

N=8 碟形算法圖

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

每一個(gè) stage 的計(jì)算順序：

• 取輸入
• 乘以轉(zhuǎn)換因子
• for section_num, for butterfly_num，執(zhí)行radixN_kernel
• 寫入輸出

看 N=8 的蝶形算法圖，stage = 1 時(shí)，運(yùn)算被分成了 4 個(gè) section，每個(gè)section的 butterfly_num = 1 。stage = 2 時(shí)， section_num = 2，butterfly_num = 2。stage = 3時(shí)， section_num = 1， butterfly_num = 4。

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

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

為當(dāng)前的 ，* sec/butterfly_stride 是每個(gè)section/butterfly* 的間隔。

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

radix-2 的冪次優(yōu)化

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

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

其中 ~ 為乘以旋轉(zhuǎn)因子 后的輸入

當(dāng) radix = 2 時(shí)，由于 , radix_2 的 DFT 矩陣形式可以寫為：

當(dāng) radix = 4 時(shí)，由于 ，radix_4的DFT 矩陣形式可以寫為：

同理推得到 radix_8 的 kernel 為：

我們先來看訪存, 現(xiàn)代處理器對于計(jì)算性能的優(yōu)化要優(yōu)于對于訪存的優(yōu)化, 在計(jì)算和訪存相近的場景下， 訪存通常是性能瓶頸。

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

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

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

以 r-4 算法為例, kernel 部分的計(jì)算為：

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

尺寸：1 * 2048（r2c) 環(huán)境：SM8150大核

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

radix-非2的冪次優(yōu)化

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

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

在復(fù)平面上根據(jù)x軸對稱，有相同的實(shí)部和相反的虛部。根據(jù)這個(gè)性質(zhì)。如下圖所示，對于每一個(gè) stage，可以合并公共項(xiàng)A，B，C，D，再根據(jù)公共項(xiàng)計(jì)算出該 stage 的輸出。

這種算法減少了很多重復(fù)的運(yùn)算。同時(shí)，在 stage>=2 的時(shí)候，同樣對 butterfly 做循環(huán)展開加并行化，進(jìn)一步減少計(jì)算的開銷。

radix-5 的優(yōu)化思想可以外推至 radix-N 。對于 radix_N 的每一個(gè) stage,計(jì)算流程為：

• 取輸入
• 乘以對應(yīng)的轉(zhuǎn)換因子
• 計(jì)算公共項(xiàng)， radix_N 有 N-1個(gè)公共項(xiàng)
• 執(zhí)行并行化的 radix_N_kernel
• 寫入輸出

其他優(yōu)化

上述兩個(gè)章節(jié)描述的是 DFT_1D 的通用優(yōu)化，在此基礎(chǔ)上還可以做更細(xì)致的優(yōu)化，可以參考本文引用的論文。

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

DFT 延展算法的原理及優(yōu)化

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

DCT

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

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

基于 DFT 的 DCT 算法流程為：

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

我們嘗試推導(dǎo)一下這個(gè)算法：

對 y[n] 依照 DFT 公式展開，整理展開的兩項(xiàng)并提取公共項(xiàng) , 根據(jù)歐拉公式和誘導(dǎo)函數(shù)，整理非公共項(xiàng) 。可以看出得到的結(jié)果正是 x[k] 和與 k 有關(guān)的系數(shù)的乘積。這樣就可以通過先計(jì)算 得到 x[n] 的 DCT 輸出 。

在理解算法的基礎(chǔ)上，我們對 DFT_1D 的優(yōu)化可以完整地應(yīng)用到 DCT 上。DCT_2D 的計(jì)算過程是依次對行、列做 DCT_1D, 我們用多線程對 DCT_1D 進(jìn)行并行，可以進(jìn)一步優(yōu)化算法。

FFT_conv

Conv 是深度學(xué)習(xí)最常見的運(yùn)算，計(jì)算conv常用的方法有 IMG2COL+GEMM， Winograd， FFT_conv。三種算法都有各自的使用場景。

FFT_conv 的數(shù)學(xué)原理是時(shí)域中的循環(huán)卷積對應(yīng)于其離散傅里葉變換的乘積. 如下圖所示， f 和 g 的卷積等同于將 f 和 g 各自做傅立葉變幻 F，進(jìn)行點(diǎn)乘并通過傅立葉逆變換計(jì)算后的結(jié)果。

直觀的理論證明可下圖。

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

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

FFT_conv 算法的流程：

• 將 Feature Map 和 Kernel 都 zero-pad 到同一個(gè)尺寸,進(jìn)行 DFT 轉(zhuǎn)換。
• 矩陣點(diǎn)乘
• 將計(jì)算結(jié)果通過 IDFT 計(jì)算出結(jié)果。

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

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

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

同時(shí)我們可以觀察到，F(xiàn)FT_conv 的核心計(jì)算模塊還是針對小圖的 DFT 運(yùn)算， 因此我們可以將前一章節(jié)對 DFT 的優(yōu)化代入此處，輔以多線程，進(jìn)一步提升 FFT_Conv 的計(jì)算效率。