注：本文為《The Scientist and Engineer‘s Guide to Digital Signal Processing》的讀書筆記系列之第八章。
Chapter 8. The Discrete Fourier Transform
- 1. The Family of Fourier Transform
- 2. Notation and Format of the Real DFT
- 3. The Frequency Domain's Independent Variable
- 4. DFT Basis Functions
- 5. Synthesis, Calculating the Inverse DFT
本章主要講述實數表示的 DFT
1. The Family of Fourier Transform
Jean Baptiste Joseph Fourier對熱的傳播很感興趣，他發表的一篇論文中提出了一個有争議的觀點：任何連續周期信号可以表示為一系列正弦信号的和。
圖1展示了如何将一個信号分解成正弦信号和餘弦信号。
圖 1-a 時域信号
圖 1-b 1-a 被分解後得到的餘弦和正弦信号
圖1-a展示了16點的時域信号，圖1-b展示了對該時域信号進行傅立葉分解得到的9個餘弦信号和9個正弦信号，每一個都具有不同的頻率和幅度。
為什麼使用正弦波而不是方波或者三角波？其實，有很多種方式都能對信号進行分解，分解的目标是得到比原始信号更容易處理的東西。正弦波和餘弦波比原始信号更簡單，因為它們具有原始信号不具有的特性：sinusoidal fidelity， 如第五章所讨論的，系統對一個正弦波的響應還是一個正弦波。隻有信号的幅度和相位可以改變，頻率和波形必須保持一緻。正弦波是唯一具有此屬性的波形。
通用術語Fourier transform可以被細分為四種類型，這是由于有四種不同類型的信号。信号可以被分為連續的和離散的，也可以被分為周期的和非周期的，總共會産生四種不同類型的信号，如圖2所示。
圖 2 四種不同類型的信号及其對應的變換名稱
連續非周期
這種信号的Fourier transform被簡單地稱為the Fourier transform。
連續周期
這種信号的Fourier transform被稱為the Fourier Series。
離散非周期
這種信号的Fourier transform被稱為the Discrete Time Fourier transform。
離散周期
這種信号的Fourier transform有時被稱為the Discrete Fourier Series，但最常被稱為the Discrete Fourier transform。
這四種類型的信号均被設定是從負無窮延伸到正無窮，即無限長。假如有一個1024點的有限長信号，那麼是不是就沒有适合其Fourier transform的版本？正弦波和餘弦波都是無限長的信号，不能使用一組無限長的信号合成有限長的信号。解決這個難題的方法是讓有限長信号look like無限長信号。這是通過想象在有限長信号的左右兩側具有無限采樣點來完成的。如果這些想象出來的無限采樣點全部是零，這時候信号看起來是離散非周期的，可以對其應用the Discrete Time Fourier transform；也可以将想象出來的采樣點認為是真實的1024個采樣點的拷貝。這時候信号看起來是離散周期的，可以對其應用the Discrete Fourier transform。
事實證明，合成一個非周期的信号需要無數個正弦波。計算機算法是不可能計算the Discrete Time Fourier transform的。可以在DSP中使用的唯一的Fourier transform類型是DFT。也就是說，數字計算機隻能使用離散且長度有限的資訊。當在理論推導或者家庭作業中思考一些數學問題時，你可能會發現你使用的是Fourier transform家族中的前三種類型，而當在使用計算機時，将僅使用DFT，本章主要講解該中類型的變換。
回顧圖1的例子，表面上看，似乎是一個16點的信号被分解成為18個正弦波，每個正弦波包含16個點。用更正式的術語說，圖1-a中所示的16點的信号必須被認為是無限長周期信号的某一周期。同樣的，圖1-b中所示的18個正弦信号中的每個都代表無限長正弦曲線中的16個點。如果将其視為由16個點正弦波合成的16個點的信号，或者視為由無限長的正弦波合成的無限長的周期信号，這真的有關系嗎？答案是：usually no, but sometimes, yes。
數學術語transform經常在數字信号進行中出現，那麼什麼是transform呢？回憶一下function，function是指将一個數或者一些數變換為另一個數的算法或者過程。transform是對function的擴充，它允許輸入輸出都是一些數。a transform is any fixed procedure that changes one chunk of data into another chunk of data。
2. Notation and Format of the Real DFT
如圖3所示，離散傅立葉變換将一個 N N N點的輸入信号變換成兩個 N / 2 + 1 N/2+1 N/2+1點的輸出信号。這兩個輸出信号包含了正弦波和餘弦波的幅度。輸入的信号一般稱為時域，輸出的信号一般稱為頻域。術語頻域被用來描述正弦波和餘弦波的幅度。
圖 3 DTF 示例
頻域包含了和時域完全一樣的資訊，隻不過是用不同的形式表示出來。标準的DSP表示法使用小寫字母表示時域信号，如： x [ ] , y [ ] x[],y[] x[],y[]和 z [ ] z[] z[]。其對應的頻域使用大寫字母表示，如 X [ ] , Y [ ] X[],Y[] X[],Y[]和 Z [ ] Z[] Z[]。為了說明，假設 N N N點的時域信号表示為 x [ ] x[] x[]，其頻域表示為 X [ ] X[] X[]，包含了實部和虛部兩部分，每部分有 N / 2 + 1 N/2+1 N/2+1個點。 X [ ] X[] X[]的實部寫作 R e X [ ] Re X[] ReX[]，虛部寫作 I m X [ ] Im X[] ImX[]，其實部包含了餘弦波的幅度，虛部包含了正弦波的幅度。
3. The Frequency Domain’s Independent Variable
圖4示範了 N = 128 N=128 N=128點的DFT。
圖 4-a 時域信号
圖 4-b、 圖 4-c 頻域信号
頻域的橫軸可以有四種不同的解釋。第一種，橫軸被标記為從0到64，使用字母 k ∈ [ 0 , N / 2 ] k\in[0, N/2] k∈[0,N/2]表示該索引，程式員最喜歡使用這種方式。圖4-b采用的這種方式。第二種，圖4-c所采用的方式，橫軸被标記為fraction of the sampling rate，這表示橫軸的取值從0到0.5，因為離散資料隻能包含介于DC和一半采樣率之間的頻率。使用字母 f f f表示該索引，且 f = ( k / N ) ∈ [ 0 , 0.5 ] f=(k/N)\in[0, 0.5] f=(k/N)∈[0,0.5]。第三種，和第二種類似，隻不過橫軸上乘了 2 π 2\pi 2π。使用字母 ω = ( 2 π f ) ∈ [ 0 , π ] \omega=(2\pi f)\in[0, \pi] ω=(2πf)∈[0,π]表示該索引，其表示natural frequency，機關為弧度。第四種，使用模拟頻率Hz表示。
4. DFT Basis Functions
DFT中所使用的正弦波和餘弦波通常稱為DFT的基函數。換句話說，DFT的輸出是一組表示幅度的數字。基函數是一組具有機關幅度的餘弦波和正弦波。DFT的基函數如下：
c k [ i ] = c o s ( 2 π k i / N ) c_k[i]=cos(2\pi ki/N) ck[i]=cos(2πki/N)
s k [ i ] = s i n ( 2 π k i / N ) (1) s_k[i]=sin(2\pi ki/N)\tag{1} sk[i]=sin(2πki/N)(1)
其中， c k [ ] c_k[] ck[]是餘弦波，其幅度存儲在 R e X [ k ] ReX[k] ReX[k]中， s k [ ] s_k[] sk[]是正弦波，其幅度存儲在 I m X [ k ] ImX[k] ImX[k]中。 k k k決定了其頻率。圖5展示了 N = 32 N=32 N=32點的DFT使用17個正弦波和17個餘弦波作為其基函數。
圖 5 DFT 的基函數
由于這些信号加起來共同構成了輸入信号，是以它們必須具有和輸入信号同樣的長度。 k k k決定了每個信号的頻率。比如， c 1 [ ] c_1[] c1[]是在N點上完成一個周期的餘弦波， c 5 [ ] c_5[] c5[]是在 N N N點上完成五個周期的餘弦波等等。這是了解基函數的一個重要概念，參數 k k k等于在 N N N點上完成的整周期數。
圖5-a展示了頻率為0的餘弦波 c 0 [ ] c_0[] c0[]，其值全為1。這表示 R e X [ 0 ] ReX[0] ReX[0]存儲時域信号所有點的平均值，也稱 R e X [ 0 ] ReX[0] ReX[0]存儲DC offset(直流偏移)。圖5-b展示了頻率為0的正弦波 s 0 [ ] s_0[] s0[]，其值全為0。 由于其不會影響時域信号的合成，是以 I m X [ 0 ] ImX [0] ImX[0]的值無關緊要，并且始終設定為零。
圖5-c和圖5-d分别展示了 c 2 [ ] c_2[] c2[]和 s 2 [ ] s_2[] s2[]，它兩均在 N N N點上完成了兩個整周期，分别對應于 R e X [ 2 ] ReX[2] ReX[2]和 I m X [ 2 ] ImX[2] ImX[2]。圖5-e和圖5-f與之類似，不過要是圖中的離散點之間沒有曲線連結，肉眼将很難看出它們是正弦波和餘弦波。
圖5-g和圖5-h分别展示了最高頻率的基函數 c 16 [ ] c_{16}[] c16[]和 s 16 [ ] s_{16}[] s16[]。離散餘弦波的值在1和-1兩者之間交替，這可以解釋為在峰值處采樣連續餘弦波。相比之下，離散正弦波的值全為0，這是由于在連續正弦波的零點處進行采樣。
有個問題：如果對 N N N個點進行DFT，而輸出了 N + 2 N+2 N+2個點，那麼額外的資訊從何而來？ 答案是：輸出中的兩個點不包含任何資訊，進而使其他 N N N個點完全獨立。不攜帶任何資訊的點是 I m X [ 0 ] ImX [0] ImX[0]和 I m X [ N / 2 ] ImX [N / 2] ImX[N/2]，其值始終為零。
5. Synthesis, Calculating the Inverse DFT
可以将合成方程寫為：
x [ i ] = ∑ k = 0 N / 2 R e X ‾ [ k ] c o s ( 2 π k i / N ) + ∑ k = 0 N / 2 I m X ‾ [ k ] s i n ( 2 π k i / N ) (2) x[i]=\sum_{k=0}^{N/2}Re\overline{X}[k]cos(2\pi ki/N)+\sum_{k=0}^{N/2}Im\overline{X}[k]sin(2\pi ki/N)\tag{2} x[i]=k=0∑N/2ReX[k]cos(2πki/N)+k=0∑N/2ImX[k]sin(2πki/N)(2)
R e X ‾ [ k ] = R e X [ k ] N / 2 Re\overline{X}[k]=\frac{ReX[k]}{N/2} ReX[k]=N/2ReX[k]
I m X ‾ [ k ] = − I m X [ k ] N / 2 Im\overline{X}[k]=-\frac{ImX[k]}{N/2} ImX[k]=−N/2ImX[k]
e x p e c t f o r t w o s p e c i a l c a s e s : expect \, for \, two \, special \, cases: expectfortwospecialcases:
R e X ‾ [ 0 ] = R e X [ 0 ] N Re\overline{X}[0]=\frac{ReX[0]}{N} ReX[0]=NReX[0]
R e X ‾ [ N / 2 ] = R e X [ N / 2 ] N (3) Re\overline{X}[N/2]=\frac{ReX[N/2]}{N}\tag{3} ReX[N/2]=NReX[N/2](3)
其中， N N N是時域信号的采樣點數， x [ i ] x[i] x[i]是合成的信号，且 i ∈ [ 0 , N − 1 ] i\in [0,N-1] i∈[0,N−1]。 R e X ‾ [ k ] Re\overline{X}[k] ReX[k]和 I m X ‾ [ k ] Im\overline{X}[k] ImX[k]分别表示合成時域信号所需的餘弦波和正弦波的幅度， R e X [ k ] Re{X}[k] ReX[k]和 I m X [ k ] Im{X}[k] ImX[k]分别表頻域的實部和虛部，記住它們的關系！！！ k ∈ [ 0 , N / 2 ] k\in [0,N/2] k∈[0,N/2]。式 2 2 2和 3 3 3一起定義了逆DFT。
表1展示了兩種實作逆DFT的編碼方法。
100 /*THE INVERSE DISCRETE FOURIER TRANSFORM
110 The time domain signal, held in XX[], is calculated from the frequency domain signals,
120 held i REX[] and IMX[].*/
130
140 float XX[512]; // XX[] holds the time domain signal
150 float REX[257]; // REX[] holds the real part of the frequency domain
160 float IMX[257]; // IMX[] holds the imaginary part of the frequency domain
170
180 const float PI = 3.14159265; // Set the constant, PI
190 const int N = 512 // N is the number of points in XX[]
200
210 // Initing REX[] and IMX[]
220
230
240 // Find the cosine and sine wave amplitudes using Eq. 3
250 for (int k=0; k <257; k++){
260 REX[k] = REX[k] / (N/2);
270 IMX[k] = -IMX[k] / (N/2);
280 }
290
300 REX[0] = REX[0] / 2;
310 REX[256] = REX[256] / 2;
320
330
340 for (int i=0; i<512; i++){ // Zero XX[] so it can be used as an accumulator
350 XX[i] = 0;
360 }
370
380 /*Method 1: Loop through each frequency generating the entire length of the sine and cosine waves, and add them to the accumulator signal, XX[] */
390 for (int k=0; k<257; k++){
400 for (int i=0; i<512; i++){
410 XX[i] = XX[i] + REX[k] * cos(2*PI*k*i/N);
420 XX[i] = XX[i] + IMX[k] * sin(2*PI*k*i/N);
430 }
440 }
450
460
470
480
490
500 /*Method 2: Loop through each sample in the time domain, and sum the corresponding samples from each cosine and sine wave*/
510 for (int i=0; i<512; i++){
520 for (int k=0; k<257; k++){
530 XX[i] = XX[i] + REX[k] * cos(2*PI*k*i/N);
540 XX[i] = XX[i] + IMX[k] * sin(2*PI*k*i/N);
550 }
560 }
圖6展示了逆DFT的操作以及頻域信号和合成時域信号所需的波形的幅度差異。
圖 6 逆 DFT
圖6-a展示了一個待合成的時域信号，是一個沖擊信号，其第一個采樣值為32，其餘采樣值均為0。圖6-b展示了該信号的頻域表示，其實部是值為32的常量，其虛部均為0（這兒沒有展示出來）。正如下一章節讨論的那樣，這是一個重要的DFT變換對：時域的沖擊信号對應于頻域的常量。現在要注意的是，圖6-b是圖6-a的DFT，圖6-a是圖6-b的逆DFT。式 3 3 3将頻域信号圖6-b轉換為餘弦波的幅度，即圖6-c。可以看出，除了第一個和最後一個餘弦波的幅度為1之外，其餘餘弦波的幅度均為2。此處沒有展示正弦波的幅度，因為，它們均為0值。式 2 2 2将餘弦波的幅度圖6-b轉換為時域信号圖6-a。
上述部分描述了頻域信号和正弦信号的幅度有何不同，但并沒有解釋為什麼會不同。它們兩者之間的差異是因為頻域信号是定義在spectral density上的。圖7展示了這是如何工作的。它展示了一個32點的時域信号所對應的頻域信号的實部。spectral density describes how much signal(amplitude) is preset per unit of bandwidth。為了将正弦信号的幅度轉化為spectral density，divide each amplitude by the bandwidth represented by each amplitude。那麼，問題來了：如何确定頻域每個離散頻率所表示的帶寬？
持續更新中。。。。。。