從傅里葉級數（Fourier series）到離散傅里葉變換（Discrete Fourier transform）

從傅里葉級數（Fourier series）到離散傅里葉變換（Discrete Fourier transform）

一. 傅里葉級數（FS）

首先從最直觀的開始，我們有一個信號\(x(t)\)（滿足Dirichelet條件），先假設它是周期的，為了研究它，我們使用級數將之展開，展開方法如下

\[x(t)=\sum_{k=0}^{\infty}a_ke^{jkw_0t}\tag{1} \]

現在問題就是如何求解\(a_k\)。因為三角函數是正交系，即

\[\forall \theta_1 \ne \theta_2 ,都有\int_{2\pi}e^{j(\theta_1 - \theta_2)t}\mathrm{d}t=0 \tag{2} \]

這樣我們就可以進行如下操作：

對（1）式左右同時積分：

\[\int_{T}x(t)e^{-jk\frac{2\pi}{T}t} \mathrm{d}t=\sum_{k'=o}^{k=\infty}\int_{T}a_ke^{j(k'-k)\frac{2\pi}{T_0}}\mathrm{d}t \tag{3} \]

由於

\[\sum_{k=o}^{k=\infty}\int_{T}a_ke^{j(k'-k)\frac{2\pi}{T_0}}\mathrm{d}t=\begin{cases} 0,& k\ne k'\\ T,&k=k' \end{cases}\tag{4} \]

故（3）積分的結果如下：

\[\int_{T}x(t)e^{-jk\frac{2\pi}{T}t} \mathrm{d}t =a_kT \]

於是我們得到了\(a_k\)：

\[a_k=\frac{1}{T}\int_{T}x(t)e^{-jk\frac{2\pi}{T}t} \mathrm{d}t \tag{5} \]

這樣我們就把一個周期信號分解成了一系列模為\(a_k\)的復周期信號的組合，形象表示就是這樣：

圖片來源：傅里葉分析之掐死傅里葉教程https://zhuanlan.zhihu.com/p/19763358（推薦大家看看）

OK，到這里我們成功將一個周期信號分解成了傅里葉級數。但是這里也有兩個問題依然困擾着我：

1.首先（相信大家也注意到了我第一句標黃的內容），非周期有限信號能不能展開成傅里葉級數呢？

比如，如果有一個函數，它是這樣的：

\[s(t) = \begin{cases} x(t), & 0\le t \le T \\ 0, & others \\ \end{cases}\\ \tag{6} x(t)是一個周期為T，滿足Dirichelet條件的函數 \]

同樣可以在0-T上做積分，並且得到\(s(t)\)的\(a_k\)和\(x(t)\)的\(a_k\)是相同的，也就是說，可以把s(t)寫成：

\[s(t)=\begin{cases} \sum_{k=0}^{\infty}a_ke^{jkw_0t},& 0\le t\le T\\ 0,& others \end{cases} \tag{7} \]

但是根據定義這樣的一個式子不叫做級數，充其量叫級數的一部分截斷。不過我們主要研究信號，就不關心這些數學概念了。那么對於這樣的一個信號，傅里葉級數還能不能表達它的頻譜特性呢？它的頻譜跟\(x(t)\)又有什么聯系呢？,這兩個問題我們在第二部分解答。

2. 三角函數系的正交性強調的是一個周期內的積分，但是沒有說一定是最小正周期，如果我們在\(2T\)、\(3T\)或者\(nT\)上積分，會發生什么呢？

回答這個問題之前，我們先對傅里葉級數做一個比較統籌的顯示，我們不妨以\(w\)作為橫坐標，以\(a_ke^{jkw_0}\)的強度（也就是\(||a_ke^{jkw_0}||=|a_k|\)）作為縱坐標，構建一個二維直角坐標系(數據大小是隨便給的)：

分析這個圖，我們知道每兩個”柱子“之間的距離應該是\(\Delta w=w_0=\frac{2\pi}{T}\)，所以，如果我們增大T，兩個數據之間的距離就會變小，如下：

可以看到，\(T\)越大,兩個數據之間的橫軸距離就越小，不妨做個猜想：當\(T\to \infty\)時，這個離散的譜就變成了連續譜。是這樣嗎？應該說對有限信號是這樣的，對無限周期信號則不是，這是為什么呢？請看下一節傅里葉變換的解讀。

二. 傅里葉變換（FT）

2.1 傅里葉變換的推導

接着前面的猜想進行分析，為了得到可靠結論，我們進行如下數學分析：

接着（5）式開展：

\[a_k=\lim_{T\to \infty}\frac{1}{T}\int_{T}x(t)e^{-jk\frac{2\pi}{T}t} \mathrm{d}t \]

用\(w_0=\frac{2\pi}{T}\)替代\(T\)，同時將\(x(t)\)表示出來:

\[\begin{split} x(t)&=\sum_{k=0}^{\infty}[\lim_{T\to \infty}\frac{1}{T}\int_{-\infty}^{+\infty}x(t)e^{-jk\frac{2\pi}{T}t} \mathrm{d}t]e^{jk\frac{2\pi}{T}t} \\ &=\sum_{k=0}^{\infty}\lim_{w_0\to 0}[\frac{w_0}{2\pi}\int_{-\infty}^{+\infty}x(t)e^{-jkw_0t} \mathrm{d}t]e^{jkw_0t}\\ \end{split} \]

當\(w_0\to 0\)時，可用\(w_0=\mathrm{d} w\)表示它，此外，對於\(kw_0\)的累加就變成了積分（累加步長\(kw_0\)無限小）,故而上式變成了：

\[x(t)=\frac{1}{2\pi}\int_{-\infty}^{+\infty}[\int_{-\infty}^{+\infty}x(t)e^{-jwt} \mathrm{d}t]e^{jwt}\mathrm{d}w \]

這樣我們實現了將一個信號分解成連續的頻譜，而這個頻譜就是：

\[X(w)=\int_{-\infty}^{+\infty}x(t)e^{-jwt} \mathrm{d}t \\ x(t)=\frac{1}{2\pi}\int_{-\infty}^{+\infty}X(w)e^{jwt}\mathrm{d}w \tag{8} \]

（8）式就是傅里葉變換的基本形式，當然有些地方可能將系數\(\frac{1}{2\pi}\)進行了一些變換,比如：

\[X(w)=\sqrt{ \frac{1}{2\pi} }\int_{-\infty}^{+\infty}x(t)e^{-jwt} \mathrm{d}t \\ x(t)=\sqrt{ \frac{1}{2\pi} }\int_{-\infty}^{+\infty}X(w)e^{jwt}\mathrm{d}w \]

2.2 周期函數的傅里葉變換

觀察（8）這個式子，我們很快又發現了另外一個問題，積分區域是無窮，而周期函數在無窮區間積分，那必然是發散的，所以我們怎么得到周期信號的傅里葉變換呢？

是這樣處理的：

1. 首先將\(x(t)\)表示成傅里葉級數：

   \[x(t)=\sum_{k} a_ke^{jkw_ot} \]

2. 對級數做傅里葉變換，也就是：

   \[X(w)=\sum_{k} F(a_ke^{jkw_ot})=a_k\sum_{k} F(e^{jkw_0t})\tag{9} \]

​ 於是問題轉化成了如何求\(e^{jkw_0t}\)的傅里葉變換，這個就要使用傅里葉變換的基本性質了：信號乘以\(e^{jw_0t}\)相當於傅里葉變換\(X(w)\)進行頻移\(w_0\)， 這個我們也可以順便證明一下：

\[若X(w)=\int_{-\infty}^{+\infty}x(t)e^{-jwt} \mathrm{d}t \\ \begin{split} 則X'(w)&=\int_{-\infty}^{+\infty}[x(t)e^{jw_0t}]e^{-jwt} \mathrm{d}t\\ &=\int_{-\infty}^{+\infty}x(t)e^{-j(w-w_0)t} \mathrm{d}t\\ &=X(w-w_0) \end{split} \]

​ 又由於：

\[1=\frac{1}{2\pi}\int_{-\infty}^{+\infty}2\pi \delta (w)e^{jwt}\mathrm{d}w \\ 故x(t)=1的傅里葉變換為F(w)= 2\pi \delta (w) \]

​ 所以：

\[F(1\cdot e^{jkw_0t})=F_{w-w_0}(1)=2\pi \delta (w-w_0)\tag{10} \]

​ 把（10）代入（9），就得到了周期函數的傅里葉變換：

\[X(w)=a_k\sum_{k} 2\pi \delta(w-kw_0)\tag{11} \]

可以看出，周期函數的傅里葉變換是一系列沖擊串，這個倒也符合直覺：

傅里葉變換每個\(X(w)\)可以理解為信號頻率為\(w\)的分量的大小，對於周期信號，它的每一個頻率分量都是一個無限長信號，所以每一個分量的能量都是無窮的，所以在頻譜上就表示為沖擊串了。

2.3 有限長非周期信號的傅里葉變換

好了，現在我們來解答第一節中的第二個問題：有限長非周期信號的傅里葉變換和周期信號頻譜的關系。

對於有限長信號：

\[s(t) = \begin{cases} x(t), & 0\le t \le T \\ 0, & others \\ \end{cases}\\ \tag{6} x(t)是一個周期為T，滿足Dirichelet條件的函數 \]

可以換一種表達方式:

\[s(t)=x(t)[u(t)-u(t-T)] \]

要求\(F[s(t)]\)，我們首先需要以下幾個結論：

1. 時域相乘等於頻域卷積，證明:

   \[若Z(w)=X(w)*Y(w),則\\ \begin{split} z(t)&=F^{-1}(Z(w))\\ &=\int_\infty [\int_\infty X(m)Y(w-m)\mathrm{d}m]e^{jwt}\mathrm{d}t \\ &=\int_\infty X(m)[\int_{\infty}Y(w-m)e^{jwt}\mathrm{d}t]\mathrm{d}m\\ &=\int_\infty X(m)y(t)e^{jmt}\mathrm{d}m\\ &=x(t)y(t)\\ 得證 \end{split} \]

   注：如果從采用\(z(t)=x(t)y(t)出發，去證明Z(w)=X(w)*Y(w)\)會比較麻煩。

2. \(F[u(t)-u(t-T)]=\frac{1-e^{-jwT}}{jw}\)，證明:

   \[\begin{split} F[u(t)-u(t-T)]&=\int_{0}^{T}x(t)e^{-jwT}\mathrm{d}t\\ &=\frac{1-e^{-jwt}}{jw} \end{split} \]

   這里的T是信號s(t)的長度，也是x(t)的最小正周期，也就是有\(T=\frac{2\pi}{w_0}\)。
   這個函數展現的頻譜圖像如下：

   

好了，我們現在可以開始計算（6）式\(s(t)\)的傅里葉變換了：

\[F(s(t))=F(x(t))*F[u(t)-u(t-T)]\\ 在(11)有X(w)=a_k\sum_{k} 2\pi \delta(w-kw_0)\\ 故F(s(t))=2\pi a_k\sum_{k}\delta(w-kw_0)*\frac{1-e^{-jwt}}{jw}\\ 又有\delta (t-t_0)*x(t)=x(t-t_0) \]

及\(\sum_{k}\delta(w-kw_0)*\frac{1-e^{-jwT}}{jw}\)就是\(F[u(t)-u(t-T)]\)以\(kw_0\)移位疊加,下面一組圖表示這個過程：

簡單起見，我們令|X(w)|為：

對應的s(t)的頻譜，就是將下面兩個信號疊加：

於是得到|S(w)|：

可以看到：在周期信號的頻譜里的沖擊串變得平滑了，或者換句話說，原來集中在某個\(kw_0\)頻率的能量一部分泄露到了全頻域。

這個在物理意義上也比較好理解，把周期信號截斷了，實際上是把一部分信號值變成了0，從某個值突然變到0(x(T)->0)，這里沒有平滑過渡，是一個突變，而我們知道，突變包含所有頻率分量。

有意思的是，我們大多數時候采集的信號都是截斷的————截斷后是工程問題，截斷前是數學問題。至此，我們完成了傅里葉級數到傅里葉變換的推導。

三、時域離散化

在實際應用中，我們大多處理的是數字信號，也就是時域離散的信號，我們來看看時域離散對信號頻譜的影響。

這個時候，我們可以把信號表示為

\[x(n)=x(t)\sum_{n=-\infty}^{+\infty}\delta(t-nT)\\ x(t)是連續信號，T是采樣周期（即每隔T時間從x(t)上取一個點） \]

根據（10），利用傅里葉變換的對偶性，可以得出：

\[F(\sum_{n=-\infty}^{+\infty}\delta(t-nT))=w_0\sum_{k=-\infty}^{+\infty}\delta(w-kw_0)\\ w_s=\frac{2\pi}{T} \]

所以\(x(t)\)的傅里葉變換為：

\[F(x(n))=w_s\sum X(w)*\delta (w-kw_s)=w_s\sum X(w-kw_s)\\ X(w)是x(t)的傅里葉變換 \]

也就是說，\(x(n)\)的傅里葉變換就是下\(x(t)\)的傅里葉變換以\(w_s\)為周期，移位疊加。

如上圖，假設中間的三角形就是\(X(w)\)，則它經過移位疊加后得到的上圖就是\(F(x(n))\)。

從這張圖我們也可以看出一個重要的定理——奈奎斯特采樣定理：采樣頻率需要高於信號最高頻率的兩倍，因為如果不能滿足

\[w_m < w_s-w_m,即w_s>2w_m \]

就會發生頻譜混疊（上圖的兩個相鄰三角形疊在一起了）。

這樣我們就得到了時域離散信號的傅里葉變換，但是顯然，這樣得到的頻譜也是連續的，計算機存儲的是離散信號，那么如何將頻域也離散化呢？請看下節。

四、DFT

回顧前面的內容，我們發現，滿足頻域離散的只有一種信號，那就是周期信號。於是讓信號離散化的方法就呼之欲出了：

假設我們采集到的N個數據其實是一個周期信號的一個周期，或者換句話說，我們把這N點數據以N為周期進行延拓，那么接下來要做的就是對一個周期信號的傅里葉變換了，又時域是離散的，傅里葉變換就退化成了傅里葉級數形式：

\[F(x(n))=\sum_{n=0}^{N-1}x(n)W_N^{nk}\\ W_{N}^{nk}=e^{-j\frac{2\pi}{N}nk}（為什么要搞出個W_{N}^{nk}呢？我也不知道） \]

現在好了，時域和頻域都是離散的了，而且頻域和時域一樣只有N個點。並且通過上面時域離散造成頻域以采樣角頻率\(w_s\)拓展的結論，我們可以知道，這個數字頻譜，其實就是在長度為\(w_s\)的區間內插入了N個點，所以沒兩個點之間的頻率間隔為

\[\Delta w=\frac{w_s}{N} \]

這個\(\Delta w\)也叫做頻率分辨率，也就是通過這個頻譜能區分的兩個頻率最小的間隔。

這里我們再對模擬角頻率\(\Omega\)和數字角頻率做一個區分：

在模擬域，頻率\(f\)表示的是一個信號每秒重復變化的次數，模擬角頻率\(\Omega=2\pi f\)表示的就是單位時間內信號相位角變換的弧度。可以把一個重復的信號想象成一個點做圓周運動，\(f\)表示它單位時間繞了多少圈，\(w\)表示它單位時間轉過多少弧度。

在數字域，我們同樣可以通過單位時間信號相位角的該變量來表示角頻率，但是數字域只有單位1，無連續時間的概念，長度為N的序列，兩個點相隔的時間是\(\frac{1}{f_s}\)，所以數字域的“時間”每改變1，角度改變應該是：

\[\begin{split} w&=\Omega \frac{1}{f_s}\\ &=\Omega T \end{split} \]

上式還可以寫作\(w=2\pi \frac{f}{f_s}\)，又由采樣定理知\(f_s>2f\),所以\(w<\pi\)，但是我們常常取數字頻域的\(0-f_s\)進行研究，也就是說\(w\)最大取到\(2\pi\)(實際上是0移位一個周期得到的)，可見：\(w\)被限制在了[0 \(2\pi\)]，所以我們也稱數字角頻率\(w\)為歸一化頻率。

---

其實上述所有的推導都能在《信號與系統》和《數字信號處理》教材里找到，我只是做了一個整理和總結，建議大家還是精讀課本，寫書的人一般比寫博客的人強。