The Fast Fourier Transform快速傅立葉轉換高速フーリエ変換

15.1 Computing the DFT Directly15.1 直接計算 DFT15.1 DFT を直接計算する

15.1.1 The DFT and twiddle factors15.1.1 DFT 與旋轉因子15.1.1 DFT と回転因子

Chapter 8 introduced the DFT as the discrete transform of a discrete (finite-length) sequence. Recall its definition, written compactly using the twiddle factor（旋轉因子） $W_N$:第 8 章把 DFT 介紹成離散（有限長度）序列的離散轉換。回想它的定義，用旋轉因子 $W_N$ 寫成緊湊的形式：第 8 章では DFT を、離散（有限長）系列に対する離散変換として導入した。その定義を回転因子 $W_N$ を使って簡潔に書き直すと次のようになる。

$W_N$ is a unit-magnitude complex number sitting at angle $-2\pi/N$ on the unit circle. Because $W_N^{r}$ only depends on $r \bmod N$, the powers $W_N^r$ take on only $N$ distinct values — the $N$-th roots of unity. Two consequences of this structure drive everything in this chapter:$W_N$ 是單位圓上角度為 $-2\pi/N$ 的單位長度複數。由於 $W_N^{r}$ 只取決於 $r \bmod N$，冪次 $W_N^r$ 只會取到$N$ 個相異值——也就是 $N$ 次單位根。這個結構帶來的兩個推論，支配了本章的一切：$W_N$ は単位円上の角度 $-2\pi/N$ に位置する、大きさ 1 の複素数である。$W_N^{r}$ は $r \bmod N$ にしか依存しないので、べき $W_N^r$ は$N$ 個の相異なる値——すなわち 1 の $N$ 乗根——しか取らない。この構造から導かれる二つの帰結が、本章のすべてを動かしている。

In matrix form the DFT is a dense $N\times N$ matrix–vector product,寫成矩陣形式，DFT 就是一個稠密的 $N\times N$ 矩陣與向量的乘積，行列の形で書けば、DFT は密な $N\times N$ 行列とベクトルの積である。

but the matrix is anything but arbitrary: every entry is one of only $N$ values, and the entries are arranged in a highly regular pattern. Figure 15-1 visualizes the real and imaginary parts of the $N=32$ DFT matrix — the structure（規律性結構）is exactly the opportunity for efficiency that the FFT exploits.但這個矩陣絕非任意：每個元素都只是那 $N$ 個值之一，而且排列得極為規律。圖 15-1 把 $N=32$ 的 DFT 矩陣的實部與虛部畫了出來——這份結構正是 FFT 所利用的效率機會。しかしこの行列は決して任意ではない。各要素はあの $N$ 個の値のどれかにすぎず、しかも極めて規則的に並んでいる。図 15-1 は $N=32$ の DFT 行列の実部と虚部を可視化したものである。この構造こそ、FFT が利用する効率化の余地にほかならない。

15.1.2 Cost of direct evaluation15.1.2 直接計算的代價15.1.2 直接計算のコスト

Evaluate the DFT sum literally. For each output point $X[k]$ we need $N$ complex multiplies and $N-1$ complex adds; there are $N$ output points ($k = 0,\dots,N-1$). Converting to real arithmetic with照字面把 DFT 的和算出來。對每一個輸出點 $X[k]$ 都需要 $N$ 次複數乘法與 $N-1$ 次複數加法；而輸出點共有 $N$ 個（$k = 0,\dots,N-1$）。再用下式換算成實數運算DFT の和を文字どおりに計算してみる。一つの出力点 $X[k]$ ごとに複素乗算が $N$ 回、複素加算が $N-1$ 回必要であり、出力点は $N$ 個ある（$k = 0,\dots,N-1$）。次式で実数演算に換算すると

$(a+jb)(c+jd) = (ac-bd) + j(ad+bc)$ = 4 real mults + 2 real adds,   1 complex add = 2 real adds,$(a+jb)(c+jd) = (ac-bd) + j(ad+bc)$ = 4 次實數乘法 + 2 次實數加法，  1 次複數加法 = 2 次實數加法，$(a+jb)(c+jd) = (ac-bd) + j(ad+bc)$ = 実数乗算 4 回 + 実数加算 2 回、  複素加算 1 回 = 実数加算 2 回、

the total damage is:總開銷就是：総計は次のようになる。

15.1.3 Goertzel's algorithm15.1.3 Goertzel 演算法15.1.3 Goertzel のアルゴリズム

Before attacking the full DFT, the lecture gives a first taste of exploiting twiddle structure: Goertzel's algorithm（Goertzel 演算法） computes a single DFT bin as the output of a small recursive filter. Start from the DFT sum and pull out a factor $W_N^{kN} = 1$:在正面處理整個 DFT 之前，課程先讓我們嚐一口利用旋轉因子結構的滋味：Goertzel 演算法把單一個 DFT 頻格算成一個小型遞迴濾波器的輸出。從 DFT 的和出發，提出一個因子 $W_N^{kN} = 1$：完全な DFT に取りかかる前に、講義は回転因子の構造を利用する最初の一例を示す。Goertzel のアルゴリズムは、一つの DFT ビンを小さな再帰フィルタの出力として計算する。DFT の和から出発し、因子 $W_N^{kN} = 1$ をくくり出す。

The right-hand sum looks like a convolution evaluated at time $N$: it is $\sum_\ell x[\ell]\, h_k[N-\ell]$ with $h_k[n] = W_N^{-kn}u[n]$. Therefore右邊那個和看起來就像摺積在時刻 $N$ 的值：它就是 $\sum_\ell x[\ell]\, h_k[N-\ell]$，其中 $h_k[n] = W_N^{-kn}u[n]$。因此右辺の和は、時刻 $N$ で評価した畳み込みの形をしている。すなわち $h_k[n] = W_N^{-kn}u[n]$ とおいた $\sum_\ell x[\ell]\, h_k[N-\ell]$ である。したがって

So one first-order complex recursion, clocked $N+1$ times (feeding $x_e[N]=0$ at the end), delivers $y_k[N] = X[k]$:所以只要一條一階複數遞迴，跑 $N+1$ 拍（最後一拍餵入 $x_e[N]=0$），就能得到 $y_k[N] = X[k]$：したがって 1 次の複素再帰を $N+1$ 回クロックし（最後に $x_e[N]=0$ を入力する）だけで、$y_k[N] = X[k]$ が得られる。

Properties claimed on the slide, worth memorizing:投影片上列出的性質，值得記下來：スライドが挙げている性質で、覚えておく価値があるもの。

• Separate "filters" for each $X[k]$ — you can compute just a few bins of interest (e.g. DTMF tone detection) without doing a whole DFT.每個 $X[k]$ 各自一個「濾波器」——只要算你關心的少數幾個頻格（例如 DTMF 按鍵音偵測），不必做整個 DFT。$X[k]$ ごとに別々の「フィルタ」——関心のあるいくつかのビンだけを計算すればよく（たとえば DTMF トーン検出）、DFT 全体を実行する必要はない。
• No large buffer, no coefficient table — each filter needs one stored coefficient $W_N^{-k}$ and one state variable; samples are consumed as they arrive.不需要大緩衝區，也不需要係數表——每個濾波器只要存一個係數 $W_N^{-k}$ 與一個狀態變數；樣本一到就吃掉。大きなバッファも係数表も不要——各フィルタが保持するのは係数 $W_N^{-k}$ 一つと状態変数一つだけであり、標本は到着したそばから消費される。
• Same complexity as the direct full DFT if you do compute all $N$ bins: $4N^2$ mults, $4N^2 - 2N$ adds.若真的把全部 $N$ 個頻格都算出來，複雜度就和直接計算整個 DFT 一樣：$4N^2$ 次乘法、$4N^2 - 2N$ 次加法。$N$ 個のビンをすべて計算するなら、計算量は直接 DFT と同じである。乗算 $4N^2$ 回、加算 $4N^2 - 2N$ 回。
• Halving the multiplies: make the denominator real by multiplying numerator and denominator by $(1 - W_N^{k}z^{-1})$:把乘法減半：分子分母同乘 $(1 - W_N^{k}z^{-1})$，讓分母變成實數：乗算を半分にする：分子と分母に $(1 - W_N^{k}z^{-1})$ を掛けて、分母を実数にする。

The recursion now uses the real coefficient $2\cos(2\pi k/N)$ — only 2 real mults per step — and the complex numerator $(1 - W_N^k z^{-1})$ is evaluated only once, at the last step $n=N$.現在這條遞迴用的是實數係數 $2\cos(2\pi k/N)$——每一步只要2 次實數乘法——而複數的分子 $(1 - W_N^k z^{-1})$ 只在最後一步 $n=N$ 算一次。この再帰が使うのは実数係数 $2\cos(2\pi k/N)$ であり、1 ステップあたり実数乗算 2 回だけで済む。複素の分子 $(1 - W_N^k z^{-1})$ は最後のステップで一度だけ、$n=N$ において評価される。

15.2 Decimation-in-Time FFT15.2 時間抽取 FFT15.2 時間間引き FFT

15.2.1 Splitting the DFT into even/odd halves15.2.1 把 DFT 拆成偶數／奇數兩半15.2.1 DFT を偶数・奇数の半分に分割する

The decimation-in-time（時間抽取）DIT FFT rearranges the DFT sum into two halves: terms with even time index $n=2m$ and terms with odd index $n = 2m+1$ (assume $N$ even, ultimately $N = 2^M$):時間抽取（DIT）FFT 把 DFT 的和重新排成兩半：時間索引為偶數 $n=2m$ 的項，以及索引為奇數 $n = 2m+1$ 的項（先假設 $N$ 是偶數，最終取 $N = 2^M$）：時間間引き（DIT）FFT は、DFT の和を二つの半分に並べ替える。時間索引が偶数 $n=2m$ の項と、奇数 $n = 2m+1$ の項である（$N$ は偶数、最終的には $N = 2^M$ とする）。

The key step uses $W_N^{2mk} = W_{N/2}^{mk}$ (index halving). $X_0$ is the $N/2$-point DFT of the even-indexed samples, and $X_1$ is the $N/2$-point DFT of the odd-indexed samples. Because an $N/2$-point DFT is periodic in $k$ with period $N/2$, we evaluate both at $\langle k\rangle_{N/2}$ ($k$ modulo $N/2$). The result, in the lecture's compact notation:關鍵的一步用到 $W_N^{2mk} = W_{N/2}^{mk}$（索引減半）。$X_0$ 是偶數索引樣本的 $N/2$ 點 DFT，$X_1$ 則是奇數索引樣本的 $N/2$ 點 DFT。由於 $N/2$ 點 DFT 對 $k$ 具有週期 $N/2$ 的週期性，兩者都在 $\langle k\rangle_{N/2}$（也就是 $k$ 對 $N/2$ 取模）處取值。用課程的緊湊記法寫出來，結果是：鍵となる一歩で $W_N^{2mk} = W_{N/2}^{mk}$（索引の半減）を使う。$X_0$ は偶数索引の標本の $N/2$ 点 DFTであり、$X_1$ は奇数索引の標本の $N/2$ 点 DFTである。$N/2$ 点 DFT は $k$ について周期 $N/2$ をもつので、どちらも $\langle k\rangle_{N/2}$（$k$ を $N/2$ で割った余り）で評価する。講義の簡潔な記法で書けば、結果は次のとおりである。

Why does this save anything? Count again:這樣到底省在哪裡？再數一次：これでなぜ節約になるのか。もう一度数えてみる。

One split halves the computation (the $\varepsilon$ accounts for the $N$ extra twiddle multiplies and adds needed to combine the halves). We can evaluate an $N$-pt DFT as two $N/2$-pt DFTs plus a few mults/adds.拆一次就把計算量減半（其中 $\varepsilon$ 代表合併兩半時額外需要的 $N$ 次旋轉因子乘法與加法）。我們可以把一個 $N$ 點 DFT 算成兩個 $N/2$ 點 DFT 再加上少數幾次乘法／加法。1 回の分割で計算量は半分になる（$\varepsilon$ は、二つの半分を合成するために余分に必要な $N$ 回の回転因子乗算と加算を表す）。$N$ 点 DFT は、二つの $N/2$ 点 DFT と少数の乗算・加算で評価できるのである。

15.2.2 One-stage DIT flowgraph15.2.2 單級 DIT 訊號流程圖15.2.2 1 段の DIT 信号流れ図

Writing the combination for all $k = 0,\dots,N-1$ and using periodicity $X_0[\langle k\rangle_{N/2}]$, $X_1[\langle k\rangle_{N/2}]$, the structure for $N=8$ is shown in Fig. 15-3. Note carefully — as the slide stresses — the twiddle factors $W_N^k$ always apply to the odd-half output $X_1[\cdot]$, and they are NOT mirror-imaged: output $X[k+N/2]$ uses $W_N^{k+N/2} = -W_N^k$, not $W_N^{k}$ again.把合併式對所有 $k = 0,\dots,N-1$ 寫出來，並利用週期性 $X_0[\langle k\rangle_{N/2}]$、$X_1[\langle k\rangle_{N/2}]$，$N=8$ 時的結構就如圖 15-3 所示。要特別注意——投影片也強調過——旋轉因子 $W_N^k$ 永遠只乘在奇數那一半的輸出 $X_1[\cdot]$ 上，而且它們並不是鏡像對稱的：輸出 $X[k+N/2]$ 用的是 $W_N^{k+N/2} = -W_N^k$，而不是再用一次 $W_N^{k}$。合成式をすべての $k = 0,\dots,N-1$ について書き下し、周期性 $X_0[\langle k\rangle_{N/2}]$、$X_1[\langle k\rangle_{N/2}]$ を使うと、$N=8$ の場合の構造は図 15-3 のようになる。スライドが強調するとおり、ここは注意を要する。回転因子 $W_N^k$ は常に奇数側の出力 $X_1[\cdot]$ にのみ掛かり、しかも鏡映にはならない。出力 $X[k+N/2]$ が使うのは $W_N^{k+N/2} = -W_N^k$ であって、$W_N^{k}$ をもう一度使うのではない。

15.2.3 Recursing to multiple stages15.2.3 遞迴到多級15.2.3 多段への再帰

If decomposing one $\mathrm{DFT}_N$ into two smaller $\mathrm{DFT}_{N/2}$'s speeds things up — why not further divide into $\mathrm{DFT}_{N/4}$'s? Apply the identical split to $X_0[k]$ (the even-sample DFT): split its inputs into even-of-even and odd-of-even subsets:如果把一個 $\mathrm{DFT}_N$ 分解成兩個較小的 $\mathrm{DFT}_{N/2}$ 可以加速——那為什麼不再往下分成 $\mathrm{DFT}_{N/4}$ 呢？把完全相同的拆分手法套用到 $X_0[k]$（偶數樣本的那個 DFT）：把它的輸入再拆成偶中之偶與偶中之奇兩個子集：一つの $\mathrm{DFT}_N$ を二つの小さな $\mathrm{DFT}_{N/2}$ に分解して速くなるなら、さらに $\mathrm{DFT}_{N/4}$ まで分割してはどうか。まったく同じ分割を $X_0[k]$（偶数標本の DFT）に適用し、その入力を偶数の中の偶数と偶数の中の奇数の部分集合に分ける。

where $X_{00}$ is the $N/4$-pt DFT of the even points within the even subset of $x[n]$ (i.e. $x[0],x[4],x[8],\dots$), $X_{01}$ the $N/4$-pt DFT of the odd points of the even subset, and so on. Note that the second-level twiddle is $W_{N/2}^{k}$ — for $N=8$ that means $W_4^k = W_8^{2k}$, so only $W_8^0$ and $W_8^2$ appear in the second stage.其中 $X_{00}$ 是 $x[n]$ 偶數子集之內那些偶數點的 $N/4$ 點 DFT（也就是 $x[0],x[4],x[8],\dots$），$X_{01}$ 則是偶數子集中奇數點的 $N/4$ 點 DFT，其餘依此類推。注意第二層的旋轉因子是 $W_{N/2}^{k}$——$N=8$ 時就是 $W_4^k = W_8^{2k}$，所以第二級只會出現 $W_8^0$ 與 $W_8^2$。ここで $X_{00}$ は $x[n]$ の偶数部分集合の中にある偶数点の $N/4$ 点 DFT（すなわち $x[0],x[4],x[8],\dots$）であり、$X_{01}$ は偶数部分集合の中の奇数点の $N/4$ 点 DFT、以下同様である。2 段目の回転因子が $W_{N/2}^{k}$ である点に注意せよ。$N=8$ ならこれは $W_4^k = W_8^{2k}$ であり、2 段目には $W_8^0$ と $W_8^2$ しか現れない。

Keep recursing until we hit 2-point DFTs, which need no multiplies at all:一直遞迴下去，直到碰上2 點 DFT，它根本不需要任何乘法：そのまま再帰を続け、2 点 DFT に行き着く。これは乗算をまったく必要としない。

An $N = 2^M$-point DFT thus reduces to $M = \log_2 N$ stages of twiddle factors and summation — the $O(N^2)$ part vanishes entirely.於是一個 $N = 2^M$ 點 DFT 就化約成 $M = \log_2 N$ 級的旋轉因子與加總——$O(N^2)$ 那一部分完全消失了。こうして $N = 2^M$ 点 DFT は、回転因子と総和からなる $M = \log_2 N$ 段へと帰着し、$O(N^2)$ の部分は完全に消え去る。

15.2.4 The butterfly and its simplification15.2.4 蝶形運算與它的化簡15.2.4 バタフライ演算とその簡略化

At every stage the elementary operation pairs two values $r$ and $r + N/2$ apart and produces two outputs. This is the butterfly（蝶形運算）. As first derived, it needs two complex multiplies, $W_N^r$ and $W_N^{r+N/2}$. But by half-period symmetry,在每一級，基本運算都把相隔 $r$ 與 $r + N/2$ 的兩個值配成一對，產生兩個輸出。這就是蝶形運算。照最初推導出來的樣子，它需要兩次複數乘法，$W_N^r$ 與 $W_N^{r+N/2}$。但由半週期對稱性，どの段でも、基本演算は $r$ と $r + N/2$ だけ離れた二つの値を組にし、二つの出力を生む。これがバタフライ演算である。最初に導いた形では複素乗算が二つ、すなわち $W_N^r$ と $W_N^{r+N/2}$ が必要になる。しかし半周期の対称性により、

so the second multiply is just the negative of the first — i.e. subtract rather than add. Each butterfly then costs one complex multiply and two complex adds:所以第二次乘法只是第一次的相反數——也就是改成相減而不是相加。於是每個蝶形只花一次複數乘法與兩次複數加法：したがって二つ目の乗算は一つ目の符号を変えたものにすぎない。つまり加えるかわりに引けばよい。こうして各バタフライのコストは複素乗算一つと複素加算二つになる。

15.3 The Complete Radix-2 FFT15.3 完整的基 2 FFT15.3 完全な基数 2 FFT

15.3.1 The 8-point DIT flowgraph15.3.1 8 點 DIT 訊號流程圖15.3.1 8 点 DIT 信号流れ図

Assembling all $M = \log_2 8 = 3$ stages of simplified butterflies gives the classic flowgraph of Fig. 15-5. Trace any output back: stage 1 computes four 2-pt DFTs, stage 2 merges them into two 4-pt DFTs (twiddles $W_8^0, W_8^2$ — i.e. $W_4^0, W_4^1$), stage 3 merges those into the final 8-pt DFT (twiddles $W_8^0,\dots,W_8^3$).把全部 $M = \log_2 8 = 3$ 級的化簡蝶形組裝起來，就得到圖 15-5 那張經典的訊號流程圖。隨便挑一個輸出往回追：第 1 級算出四個 2 點 DFT，第 2 級把它們併成兩個 4 點 DFT（旋轉因子 $W_8^0, W_8^2$——也就是 $W_4^0, W_4^1$），第 3 級再把這兩個併成最終的 8 點 DFT（旋轉因子 $W_8^0,\dots,W_8^3$）。簡略化したバタフライを $M = \log_2 8 = 3$ 段ぶん組み上げると、図 15-5 の古典的な信号流れ図が得られる。どの出力でもよいので遡ってみよ。第 1 段は 2 点 DFT を四つ計算し、第 2 段はそれらを二つの 4 点 DFT に併合し（回転因子 $W_8^0, W_8^2$、すなわち $W_4^0, W_4^1$）、第 3 段がそれらを最終的な 8 点 DFT に併合する（回転因子 $W_8^0,\dots,W_8^3$）。

15.3.2 Bit-reversed input ordering15.3.2 輸入的位元反轉排序15.3.2 入力のビット反転順序

Each DIT split sends even-indexed samples up and odd-indexed samples down — i.e. it sorts on the least significant bit (LSB) of the time index. Recursing sorts on the next bit, and so on. After $M$ levels, the sample at input row $r$ is $x[\text{bit-reverse}(r)]$: the address with its $M$-bit binary representation reversed（位元反轉）.每一次 DIT 拆分都把偶數索引的樣本送上去、奇數索引的樣本送下來——也就是依時間索引的最低有效位元（LSB）做排序。再遞迴一層就依下一個位元排序，依此類推。經過 $M$ 層之後，輸入第 $r$ 列上的樣本就是 $x[\text{bit-reverse}(r)]$：也就是把位址的 $M$ 位元二進位表示反轉過來的那個位址。DIT の分割は毎回、偶数索引の標本を上へ、奇数索引の標本を下へ送る。つまり時間索引の最下位ビット（LSB）で並べ替えているのである。再帰すれば次のビットで並べ替え、以下同様となる。$M$ 段を経たあと、入力の第 $r$ 行にある標本は $x[\text{bit-reverse}(r)]$ である。すなわちアドレスの $M$ ビット 2 進表現を逆順にしたものである。

15.3.3 In-place computation15.3.3 原位計算15.3.3 その場計算

Look again at Fig. 15-5: every butterfly takes two values from rows $(p, q)$ and writes its two outputs back to the same two rows $(p, q)$. No other butterfly in that stage touches those rows. Consequently the whole FFT is an in-place algorithm（原位演算法）: one length-$N$ complex array suffices — after the initial bit-reversal shuffle, the $M$ stages overwrite the array sequentially, and the final array contents are $X[0],\dots,X[N-1]$ in natural order. Memory cost: $N$ complex words plus a small twiddle table (or none for Goertzel-style on-the-fly generation).再看一次圖 15-5：每個蝶形從第 $(p, q)$ 兩列取兩個值，又把兩個輸出寫回同樣那兩列 $(p, q)$。該級中沒有別的蝶形會碰到這兩列。因此整個 FFT 是一個原位演算法：一個長度為 $N$ 的複數陣列就夠了——在最初的位元反轉重排之後，$M$ 個級依序覆寫這個陣列，最後陣列的內容就是自然順序的 $X[0],\dots,X[N-1]$。記憶體代價：$N$ 個複數字組，再加上一張小小的旋轉因子表（若像 Goertzel 那樣即時產生，連表都不用）。図 15-5 をもう一度見てほしい。どのバタフライも第 $(p, q)$ 行から二つの値を取り、二つの出力を同じその二行 $(p, q)$ に書き戻す。その段のほかのバタフライがこの二行に触れることはない。したがって FFT 全体はその場アルゴリズムである。長さ $N$ の複素配列が一つあれば足りる。最初のビット反転の並べ替えのあと、$M$ 段が順に配列を上書きしていき、最終的な配列の中身が自然順序の $X[0],\dots,X[N-1]$ となる。メモリのコストは複素数 $N$ 語と、小さな回転因子表（Goertzel 流にその都度生成するなら表も不要）である。

15.3.4 Complexity: (N/2) log₂ N15.3.4 複雜度：(N/2) log₂ N15.3.4 計算量：(N/2) log₂ N

An $N = 2^M$-pt DFT reduces to $M = \log_2 N$ stages. Each stage has $N/2$ butterflies, each with at most 1 complex multiply (4 real mults) and 2 complex adds (4 real adds):一個 $N = 2^M$ 點 DFT 化約成 $M = \log_2 N$ 級。每一級有 $N/2$ 個蝶形，每個蝶形最多 1 次複數乘法（4 次實數乘法）與 2 次複數加法（4 次實數加法）：$N = 2^M$ 点 DFT は $M = \log_2 N$ 段に帰着する。各段には $N/2$ 個のバタフライがあり、それぞれ複素乗算が高々 1 回（実数乗算 4 回）、複素加算が 2 回（実数加算 4 回）である。

The quadratic term is gone. Figure 15-6 shows just how dramatic the difference is — at $N = 2^{16}$ the FFT needs about 3 orders of magnitude fewer multiplies.平方項消失了。圖 15-6 顯示這個差距有多戲劇性——在 $N = 2^{16}$ 時，FFT 所需的乘法次數少了大約三個數量級。二乗の項は消え去った。図 15-6 はその差がいかに劇的かを示している。$N = 2^{16}$ では、FFT に必要な乗算回数は3 桁ほど少ない。

15.4 FFT Variations15.4 FFT 的各種變形15.4 FFT の変種

15.4.1 FFT for other values of N15.4.1 其他 N 值的 FFT15.4.1 その他の N に対する FFT

Having $N = 2^M$ meant we could divide each stage into 2 halves — the radix-2 FFT（基數 2 FFT）. The same divide-and-conquer approach works more generally:因為有 $N = 2^M$，我們才能把每一級切成兩半——這就是基 2 FFT。同樣的分治手法可以推廣得更一般：$N = 2^M$ であったからこそ、各段を 2 等分できた——これが基数 2 FFT である。同じ分割統治の手法はもっと一般に使える。

• $N = 3^M$ → radix-3 FFT (split into thirds);$N = 3^M$ → 基 3 FFT（切成三份）；$N = 3^M$ → 基数 3 FFT（3 等分する）；
• $N = 4^M$ → radix-4 FFT — a more optimized radix-2 (fewer twiddle multiplies per output);$N = 4^M$ → 基 4 FFT——更進一步最佳化的基 2（每個輸出所需的旋轉因子乘法更少）；$N = 4^M$ → 基数 4 FFT——基数 2 をさらに最適化したもの（出力あたりの回転因子乗算が少ない）；
• etc. for any prime radix;任何質數基數也都依此類推；任意の素数の基数についても同様である；
• composite $N = a\cdot b\cdot c\cdot d$ → mixed radix（混合基數）: a different $N/r$-point split at each stage;合成數 $N = a\cdot b\cdot c\cdot d$ → 混合基數：每一級用不同的 $N/r$ 點拆分；合成数 $N = a\cdot b\cdot c\cdot d$ → 混合基数：段ごとに異なる $N/r$ 点の分割を使う；
• … or just zero-pad（補零） the sequence up to the next $N = 2^M$ (recall from Ch. 8 that zero-padding only interpolates the spectrum more finely — it changes which frequencies are sampled, not the underlying DTFT).……或者乾脆把序列補零到下一個 $N = 2^M$（回想第 8 章：補零只是把頻譜內插得更細——它改變的是取樣到哪些頻率，而不是底下的 DTFT）。……あるいは単に系列を次の $N = 2^M$ までゼロ詰めすればよい（第 8 章を思い出そう。ゼロ詰めはスペクトルをより細かく内挿するだけであり、どの周波数を標本化するかを変えるのであって、背後の DTFT を変えるわけではない）。

15.4.2 Inverse FFT via the forward FFT15.4.2 用正向 FFT 做反 FFT15.4.2 順方向 FFT による逆 FFT

Recall the IDFT:回想 IDFT：IDFT を思い出そう。

The only differences from the forward DFT are the $1/N$ scale and the sign of the exponent. Conjugating twice fixes the sign:與正向 DFT 的差別只有 $1/N$ 這個比例因子和指數的符號。取兩次共軛就能修正符號：順方向 DFT との違いは $1/N$ の係数と指数の符号だけである。共役を 2 回とれば符号は直る。

The bracketed sum is exactly a forward DFT of $x'[n] = X^*[k]\big|_{k=n}$ — a "time" sequence made from the spectrum. So the recipe is: conjugate the spectrum, run the same FFT engine, conjugate the result, divide by $N$. One FFT implementation serves both directions; in a pure real flowgraph this amounts to swapping the roles of real/imaginary parts with a couple of sign flips and $1/N$ scalings at the ports.括號裡的那個和，正好就是 $x'[n] = X^*[k]\big|_{k=n}$ 的正向 DFT——一條用頻譜做出來的「時間」序列。所以做法是：把頻譜取共軛，送進同一套 FFT 引擎，結果再取共軛，然後除以 $N$。一份 FFT 實作就能服務兩個方向；在純實數的訊號流程圖裡，這相當於把實部與虛部的角色對調，再在埠口加上幾個符號反轉與 $1/N$ 的縮放。角括弧の中の和は、$x'[n] = X^*[k]\big|_{k=n}$ の順方向 DFT にほかならない。スペクトルから作った「時間」系列である。したがって手順は、スペクトルの共役をとり、同じ FFT エンジンに通し、結果の共役をとって $N$ で割る、というだけである。FFT の実装は一つで双方向に使える。純実数の信号流れ図でいえば、実部と虚部の役割を入れ替え、端子でいくつかの符号反転と $1/N$ のスケーリングを行うことに相当する。

15.4.3 DFT of real sequences — two for the price of one15.4.3 實數序列的 DFT——一次的代價做兩條15.4.3 実数系列の DFT — 一つの値段で二つ

If $x[n]$ is pure-real, the DFT wastes multiplies: by conjugate symmetry（共軛對稱）, $X[k] = X^*[-k] = X^*[N-k]$, so half the outputs are redundant. The lecture's trick: pack two real sequences into one complex FFT. Given real $x[n]$ and $w[n]$, form $y[n] = x[n] + j\,w[n]$ and take a single $N$-pt DFT $V[k] = \mathrm{DFT}\{y[n]\} = X[k] + Y[k]$, where $Y[k] = \mathrm{DFT}\{j\,w[n]\} = jW[k]$. Then如果 $x[n]$ 是純實數，DFT 就浪費了乘法：由共軛對稱性，$X[k] = X^*[-k] = X^*[N-k]$，所以有一半的輸出是多餘的。課程的技巧是：把兩條實數序列打包進一次複數 FFT。給定實數的 $x[n]$ 與 $w[n]$，組成 $y[n] = x[n] + j\,w[n]$，然後只做一個 $N$ 點 DFT $V[k] = \mathrm{DFT}\{y[n]\} = X[k] + Y[k]$，其中 $Y[k] = \mathrm{DFT}\{j\,w[n]\} = jW[k]$。於是$x[n]$ が純実数なら、DFT は乗算を無駄にしている。共役対称性により $X[k] = X^*[-k] = X^*[N-k]$ であるから、出力の半分は冗長である。講義の技は、実数系列二つを一回の複素 FFT に詰め込むことである。実数の $x[n]$ と $w[n]$ が与えられたら $y[n] = x[n] + j\,w[n]$ を作り、$N$ 点 DFT を一回だけとって $V[k] = \mathrm{DFT}\{y[n]\} = X[k] + Y[k]$ を得る。ここで $Y[k] = \mathrm{DFT}\{j\,w[n]\} = jW[k]$ である。すると

(Here $V^*[-k] = V^*[\langle N-k\rangle_N]$, and the cancellations use the conjugate symmetry of the real-input spectra $X$ and $W$.) That is: compute the DFTs of two $N$-pt real sequences with a single $N$-pt complex DFT — Exercise 5 verifies this numerically to machine precision.（這裡 $V^*[-k] = V^*[\langle N-k\rangle_N]$，而那些抵消用到的是實數輸入頻譜 $X$ 與 $W$ 的共軛對稱性。）也就是說：用單一個 $N$ 點複數 DFT 就能算出兩條 $N$ 點實數序列的 DFT——習題 5 會用數值把這件事驗證到機器精度。（ここで $V^*[-k] = V^*[\langle N-k\rangle_N]$ であり、打ち消しには実数入力のスペクトル $X$ と $W$ の共役対称性を使っている。）すなわち、$N$ 点の実数系列二つの DFT を、$N$ 点の複素 DFT 一回で計算できるのである。演習 5 ではこれを数値的に機械精度まで確かめる。

15.5 The FFT in MATLAB15.5 MATLAB 中的 FFT15.5 MATLAB での FFT

MATLAB's fft/ifft call the FFTW library — mixed-radix $O(N\log N)$ for every length, fastest for $N = 2^M$:MATLAB 的 fft／ifft 會呼叫 FFTW 函式庫——對任何長度都是混合基數的 $O(N\log N)$，而在 $N = 2^M$ 時最快：MATLAB の fft／ifft は FFTW ライブラリを呼ぶ。どんな長さでも混合基数の $O(N\log N)$ であり、$N = 2^M$ のとき最も速い。

N  = 4096;
x  = randn(1, N);
X  = fft(x);              % N-point FFT (FFTW under the hood)
xr = ifft(X);             % inverse; max(abs(xr - x)) ~ 1e-16
X2 = fft(x, 8192);        % zero-pads x to 8192 points first

% --- timing: power-of-two vs prime length ---
N1 = 2^16;  N2 = 65537;                 % 65536 vs the prime 65537
x1 = randn(1, N1);  x2 = randn(1, N2);
t1 = timeit(@() fft(x1));               % ~ a few ms
t2 = timeit(@() fft(x2));               % noticeably slower
fprintf('N = %6d : %8.3f ms\n', N1, 1e3*t1);
fprintf('N = %6d : %8.3f ms\n', N2, 1e3*t2);

Figure 15-7 shows measured run times on this machine (numpy's FFTPACK-style engine; MATLAB/FFTW behaves the same way qualitatively): both power-of-two and prime lengths scale like $N \log N$, but the prime-length constant is ~5× worse, and the direct $O(N^2)$ matrix product blows past both almost immediately.圖 15-7 顯示在這台機器上實測的執行時間（numpy 的 FFTPACK 式引擎；MATLAB/FFTW 在定性上行為相同）：2 的冪與質數長度都按 $N \log N$ 成長，但質數長度的常數差了約 5 倍，而直接做 $O(N^2)$ 的矩陣乘積幾乎立刻就把兩者遠遠甩在後頭。図 15-7 はこの計算機で実測した実行時間である（numpy の FFTPACK 系エンジン。MATLAB/FFTW も定性的には同じ振る舞いをする）。2 のべき乗長も素数長もともに $N \log N$ でスケールするが、素数長の定数は約 5 倍悪く、直接の $O(N^2)$ 行列積はほとんど即座に両者を置き去りにして伸びていく。

15.6 Summary15.6 總結15.6 まとめ

Exercises習題演習問題

N = 64;  n = 0:N-1;
x = cos(2*pi*5*n/N) + 0.5*cos(2*pi*12*n/N);   % real sequence 1
w = sin(2*pi*9*n/N) + 0.3*randn(1,N);         % real sequence 2

y  = x + 1j*w;                                % pack into ONE complex sequence
V  = fft(y);                                  % single N-pt complex FFT
Vc = conj(V(mod(N-(0:N-1), N) + 1));          % V*[_N]  (= V*[-k])

X = (V + Vc)/2;                               % spectrum of x
W = (V - Vc)/(2j);                            % spectrum of w  (-j/2 form)

errX = max(abs(X - fft(x)))                   % ~ 3.6e-15
errW = max(abs(W - fft(w)))                   % ~ 5.3e-15

subplot(2,1,1); stem(n, abs(X)); ylabel('|X[k]|');
title('Recovered from one complex FFT');
subplot(2,1,2); stem(n, abs(W)); ylabel('|W[k]|'); xlabel('k');