如何理解離散傅里葉變換及Z變換

來自專欄小潘工程師筆記70 人贊了文章

專欄前面的幾篇文章主要介紹了連續系統的一些知識，今天呢我們簡要說一下離散系統，主要是離散傅里葉變換、Z變換以及其他的一些基礎知識。

一、什麼是採樣

我們一般對模擬信號比較熟悉，它在時域是連續的，一般長成這個樣子：

模擬信號一般可以在示波器裡面進行觀測，在很久很久以前，我們需要處理的信號只有模擬信號。但是現在不一樣了，因為我們步入了新時代——數字時代，大部分信號都變成數字式了，典型的數字信號長成這個樣子：

把模擬信號變成數字信號的過程稱之為採樣。

採樣頻率定義為： f_s=frac{1}{T_s} ， T_s 為採樣時間。

採樣是一個有規律的周期性過程，也就說，採樣會引入額外的諧波分量。舉個簡單的例子，現在有一個餘弦信號，頻率為 8Hz，表達式為：

x(t)=Acos(2picdot8cdot t)

假設我們對這個餘弦信號進行採樣，採樣頻率為 f_s=20Hz ，採樣結果如下圖，其中虛線為原始信號，菱形為採樣點的數值。

現在我們將原始曲線去掉，只留下菱形的採樣點，如下圖：

那怎樣才能復原原來的曲線呢？聰明的你可能已經想到了：用曲線擬合啊，假如所有的點都能在某個曲線上，那擬合出來的曲線很有可能就是被採樣的曲線，這樣想貌似沒什麼大問題。

然而事實卻很打臉，我們極有可能模擬出如下曲線來，也是一個餘弦，但是頻域不是原先的8Hz，也不是採樣所用的20Hz，而是12Hz，模擬後的12Hz曲線如下圖，這是為什麼？

如果我們工作足夠細緻，就會發現還可以用如下頻率的曲線擬合採樣點：28Hz，32Hz，48Hz，52Hz, ... ，也就是說採樣之後信號頻譜有很多頻率，如下圖：

那麼問題就來了，其他的頻率分量到底是哪來的？

這個問題還真不簡單，我們要做幾個鋪墊以後才能回答。

二、狄拉克梳狀函數（Dirac comb）

2.1 什麼是狄拉克(Dirac)函數

在專欄「如何理解不確定性原理——不確定or測不準」那篇文章中，我們已經認識了什麼是狄拉克函數，具體參見如下鏈接：

J Pan：如何理解不確定性原理—不確定or測不準？?

zhuanlan.zhihu.com圖標

文章中提到，狄拉克 delta 函數定義為：

delta (x)={egin{cases} +infty ,& x=0\ 0,& x eq 0end{cases}}

並且滿足

int_{-infty}^{+infty}delta(x)dx=1

從定義很容易得出這個函數有如下一些性質：

  • 對稱性質： delta(t)=-delta(t)
  • 縮放性質： delta(at)=frac{1}{a}delta(t)
  • 篩選性質： int _{-infty }^{infty }x(t)delta (t - t_0)dt =x(t_0)

2.2 什麼是狄拉克梳狀函數（Dirac comb）

前面說了，狄拉克函數具有篩選性質，但是每次只能篩選出一個點，那怎樣才能有規律的篩選出更多的採樣點呢？——狄拉克梳狀函數（Dirac comb），這個函數在電子工程領域被稱為脈衝序列（impluse train）或者是採樣函數（sampling function），有些文章也稱之為Shah function。狄拉克梳狀函數的本質就是狄拉克 delta 函數用周期為 T_s 間隔的無窮多級數（多個 delta 函數合併而成），數學表達式為：

delta_s(t)=sum_{n=-infty}^{infty}{delta(t-nT_s)}

其時域波是周期為 T_s 的單位脈衝序列，也稱之為理想抽樣函數，長相如下：

由於是周期函數，我們可以按照周期函數的傅里葉級數求其頻譜，其傅里葉級數的係數為：

c_n=frac{1}{T_s}int_{-frac{1}{2}T_s}^{frac{1}{2}T_s}delta_s(t)e^{-jfrac{2pi n}{T_s}t}dt=frac{1}{T_s} , n=0,pm1,pm2,pm3,...

所以其傅里葉級數展開為：

delta_s(t)=sum_{n=-infty}^{+infty}c_ne^{jfrac{2pi n}{T_s}t}=frac{1}{T_s}sum_{n=-infty}^{+infty}e^{jfrac{2pi n}{T_s}t}

可見，其頻率值分布在 f=nfrac{2pi}{T_s} ( n=0,pm1,pm2,pm3,... )的離散頻率點，各頻率分量的幅值均為 frac{1}{T_s} ，傅里葉級數係數如下：

其實引進了衝擊函數的概念之後，對一般的周期函數也可以求其傅里葉變換：

F(omega)=mathcal{F}[delta_s(t)]=mathcal{F}left[ sum_{n=-infty}^{+infty}{c_n e^{jfrac{2pi n}{T_s}t}} ight]

稍微變形一下；

F(omega)=sum_{n=-infty}^{+infty}c_ncdotmathcal{F}left[ { e^{jfrac{2pi n}{T_s}t}} ight]=frac{2pi}{T_s} sum_{n=-infty}^{+infty}delta(omega-nfrac{2pi}{T_s})

其中 c_n 為周期函數的傅里葉級數，其頻譜如下：

可見狄拉克梳狀函數的傅里葉變換仍然是頻域的狄拉克梳狀函數，頻譜衝激串的周期是基頻是採樣頻率f_s ，衝激強度為 2pi f_s 即frac{2pi}{T_s}。用狄拉克梳狀函數表示離散信號和離散頻譜，將給信號分析和處理帶來極大的方便。

三、什麼是離散時間傅里葉變換（DTFT）

在本專欄前面的文章中，我們介紹了什麼是傅里葉變換，已經忘記傅里葉變換是怎麼回事的童鞋可複習一下這篇文章：

J Pan：傅里葉變換後面的到底有什麼小秘密?

zhuanlan.zhihu.com圖標

傅里葉變換是用來處理連續系統的，那離散系統如何處理呢？假如對 x(t) 為連續系統信號，要轉換成離散系統，就要先進行採樣，採樣頻率為 f_s ，衝擊採樣序列為：

delta_s(t)=sum_{n=-infty}^{infty}{delta(t-nT_s)}

則取樣之後的的信號為

x_s(t)=sum_{n=-infty}^{infty}x(t)delta(t-nT_s)

連續信號的傅里葉變換的定義為：

X(jomega)=int_{-infty}^{+infty}x(t)e^{-iomega t}dt

則採樣後的信號傅里葉變換為：

X_s(jomega)=int_{-infty}^{infty}left( sum_{n=-infty}^{infty}x(t)delta(t-nT_s) ight)e^{-jomega t}dt

交換一下積分與求和順序：

X_s(jomega)=sum_{n=-infty}^{infty}int_{-infty}^{infty}delta(t-nT_s)left( x(t) e^{-jomega t} ight)dt

由 delta(t) 的篩選特性 int_{-infty}^{infty}delta(x-x_0)f(x)dx=f(x_0) 可知：

X_s(jomega)=sum_{n=-infty}^{infty}x(nT_s)e^{-jomega cdot nT_s}

對數字信號而言，我們只有 x(t) 採樣後的信號 {x[n]} ，第n個採樣發生在時間 t=nT_s 。因此可以將連續信號的的傅里葉變換寫成如下形式：

X_s(jomega)=sum_{n=-infty}^{infty}{x[n]cdot e^{-jomega cdot nT_s}}

上式稱之為離散時間傅里葉變換（Discrete Time Fourier Transform），簡稱為DTFT，也就是無限長離散信號如何進行傅里葉變換。但是呢，DTFT的結果頻率 omega 還是連續的，也是不能直接應用在數字系統中。那實際是怎麼辦呢？

四、什麼是離散傅里葉變換（DFT）

我們目前處於一個數字化的世界，如果一個信號的頻譜是連續的（DTFT），那麼我們仍然沒法在機器中處理，那怎麼辦呢？——頻率也要離散化才行，所需要的技術就是離散傅里葉變換DFT（Discrete Time Fourier Transform），即具有周期特性離散信號的傅里葉級數（就是將無限長的離散限號進行截短至N個採樣點，然後將這個N個採樣點進行周期延拓 ），DFT表達式如下：

X[k]=frac{1}{N}sum_{n=0}^{N-1}{x[n]e^{-jfrac{2pi}{N}nk}}, 0leq k<N-1

x[n]=sum_{n=0}^{N-1}{X[k]e^{jfrac{2pi}{N}nk}} , 0leq n<N-1

一般的教材，都是直接給出DFT計算公式，對DFT與DTFT以及傅里葉變和傅里葉級數的關係並沒有明確，導致很多童鞋不明覺厲，總覺得心裡不是很踏實。

稍微負責任一點的書會告訴大家： left{ e^{jfrac{2pi}{N}kn} ight} ,其中0leq k <N-1 ，在 0leq n <N-1 區間內是完備正交的，然後通過這個完備正交基獲得DFT的表達式，這多少給人馬後炮的感覺，是逆向推導而不是正向設計。

今天我們提供一個思路，來看看第一個人或許受到了怎麼的啟發，發現了DFT的演算法。我們知道，周期函數可以用傅里葉級數展開，在引入了 delta 函數後，更一般的周期函數也能被展開了，根據連續周期信號的傅里葉級數：

c_n=frac{1}{T}int_{0}^{T}f(t)e^{-jfrac{2pi}{T}nt}dt

對於連續信號 x(t) 按照採樣時間 T_s 進行抽樣 N 次 delta(t-nT_s) ，並將這 N 個數值進行周期延拓，可以得到周期的離散信號，在一個周期內，其表達式為：

x_s(t)=x(t)sum_{0}^{N-1}{delta(t-nT_s)}

可以獲得離散周期信號的傅里葉級數為： X(jkomega)=frac{1}{T}int_{-T/2}^{T/2}x(t)sum_{n=0}^{N-1}{delta(t-nT_S)e^{-jfrac{2pi}{T}kt}}dt

將積分與求和調換一下順序：

X(jkomega)=frac{1}{T}sum_{n=0}^{N-1}int_{-T/2}^{T/2}x(t){delta(t-nT_S)e^{-jfrac{2pi}{T}kt}}dt

由於 delta 衝擊函數的篩選性質，上式很容易進行離散化，很顯然， T ightarrow N ，積分號內只有採樣點 sum_{0}^{N-1}nT_s 處的值被保留，所示上式可以簡化為： X[jkomega]=frac{1}{N}sum_{n=0}^{N-1}x(nT_s){e^{-jfrac{2pi}{NT_s}knTs}}=frac{1}{N}sum_{n=0}^{N-1}x[n]{e^{-jfrac{2pi}{N}kn}}

這就是離散周期信號的傅里葉變換，理論上離散周期信號的頻譜是有無窮多的，但是由於 left{ e^{jfrac{2pi}{N}kn} ight} 的周期性，一般我們只取主值區間0leq k<N-1 進行研究。

上面我們計算了離散周期信號得頻譜，即不同頻率分量下對應的係數，下面將得出一個重要的結論，很是很多童鞋迷惑的一點，那就是離散傅里葉變換隻有離散的數值，那頻率去哪了？

komega=kfrac{2pi}{T}=kfrac{2pi}{NT_s}=2pifrac{k}{N}frac{1}{T_s}=2pi f_sfrac{k}{N}

也就說：

f_k=frac{k}{N}f_s

離散周期傅里葉變換後的第 k 個數對應的頻率是 frac{k}{N}f_s ，這個結論非常重要，在後面的Z變換的性質中會用到。

說到周期離散信號的傅里葉變換DFT（discrete Fourier transform），就不得不提FFT（fast Fourier transform ），即快速傅里葉變換，畢竟FFT在工程中更常見，那它們有什麼關係呢？——這倆貨其實呢是一回事，FFT本質也是DFT，只不過是利用DFT內部蘊藏的規律的一種簡化演算法罷了。

五、什麼是Z變換

在第三部分，我們已經介紹到，DTFT的公式是

X(jomega)=sum_{n=-infty }^{infty }{x[n]e^{-jomega nT_s} }

同樣的，DTFT需要滿足絕對可和的條件，即sum_{n=-infty }^{infty }{left| x[n] ight| } <infty 。然而事實上並不是所有的離散數列都滿足這個絕對可積條件，那怎麼辦呢？在文章「從另一個角度看拉普拉斯變換」中，我們也碰到過類似的問題，信號 f(t) 要保證能進行傅里葉變換，也要滿足絕對可積條件： int_{0}^{+infty}left| f(t) ight|<+infty

我們是怎麼做的呢？——軟的不行就來硬的唄，給 f(t) 強制乘以一個衰減函數 e^{-sigma t} ，保證

int_{0}^{+infty}left| f(t)e^{-sigma t} ight|<+infty

具體請參照

J Pan：從另一個角度看拉普拉斯變換?

zhuanlan.zhihu.com圖標

對於離散數列，是不是也可以採取這個思路呢？我們不妨試一下：我們給 x[n] 乘一個指數函數e^{-sigmacdot nT_s } ，對應連續信號的衰減函數 e^{-sigma t} ，即用 nT_s 對應時間 t 。

則函數x[n]e^{-sigma nT_s} 的DTFT為：

X(jomega)=sum_{n=-infty }^{infty }{x[n]e^{-sigma nT_s} e^{-jomega nT_s} } ，

進一步化簡：

X(jomega)=sum_{n=-infty }^{infty }{x[n]e^{-(sigma+jomega) nT_s} }

這看起來太複雜了，既不符合工程師的調性，也不符合數學家的調性。我們要簡潔，簡潔，簡潔！如果我們令

z=e^{(sigma+jomega)T_s}

則DFT的表達式變為：

X(jomega)=sum_{n=-infty}^{infty}x[n]z^{-n}

是不是簡潔多了！是不是神清氣爽了！——這就是Z變換了。

如果大家記性好的話，肯定還記得我們在推導Z變換的過程中借鑒了拉普拉斯變換的很多處理方式，那它們之間是不是有某種聯繫呢？

我們知道：

s=sigma+jomega

我們又知道;

z=e^{(sigma+jomega)T_s}

顯然能得出以下結論：

z=e^{sT_s}

也就是說， z 域是 s 域的進一步映射，z 變換和拉普拉斯變換通過某種映射連接在一起。同時，由 z=e^{sT_s} 還可以看出，穩態時， s=jomega ，則 z=e^{jomega T_s} ，這個函數大家可能很熟悉了，這尼瑪就是一個單位圓啊！如果看不出來也沒關係，可以參看小潘的另外一篇文章：J Pan：被眾人膜拜的歐拉恆等式是個什麼東東？?

zhuanlan.zhihu.com圖標

我們再深入研究一下z和 s 之間的映射關係：

由 s=sigma+jomega ，當 sigma=0 時，s=jomega ，對應的是 s 域的虛軸，而此時 z=e^{jomega T} 對應的是單位圓，也就是說z 變換將 s 域的虛軸映射成 z 域的單位圓。

當 sigma>0 時，s=sigma +jomega ，對應的是 s 域的正半軸，而此時 z=e^{sigma}e^{jomega T} ，由於 e^{sigma}>1 ，也就是說此時z 變換將s 域正半軸映射到了z 域的單位圓外部。

當 sigma<0 時，s=sigma +jomega ，對應的是 s 域的負半軸，而此時 z=e^{sigma}e^{jomega T} ，由於 e^{sigma}<1 ，也就是說此時z 變換將s 域負半軸映射到了z 域的單位圓內部。

這樣做有什麼好處呢？

如果 x(t) 和 X(s) 是拉普拉斯變換對，則 x(t-T_s) 和 e^{-sT_s}X(s) 也是變換對，如果 x(t) 延時 n imes T_s ,那麼 x(t-nT_s) 的傅里葉變換對為 (e^{-sT_s})^nX(s) ，可見：在時域進行延時操作，在頻域就相當於旋轉操作， e^{-sT_s} 大量出現在具有延時特性系統的拉普拉斯變換裡面，如果用 z=e^{sT_s} 來表示的話，表達和計算都會方便很多。簡單點來說：z 就相當於定義了一個旋轉操作符，這對於具有延時採樣的離散系統特別有用。

六、什麼是採樣定理

好了，能堅持到這的童鞋都是好奇寶寶——在等我們第一部分問題的答案呢！到底其他頻率的分量是哪來的？

我們知道，兩個信號在時域相乘，在頻域相當於卷積；在時域卷積，在頻響相當於相乘，這就是卷積定理，具體參見文章

J Pan：「卷積」其實沒那麼難以理解?

zhuanlan.zhihu.com圖標

離散信號，其實就是連續信號 x(t) 與狄拉克梳狀函數（也就是採樣函數）的相乘，這就是採樣。這是時域行為，那在頻域會怎麼樣？——卷積唄！

在本文的第二部分中，我們又介紹了，狄拉克梳狀函數無論在時域還是在頻域，其形貌都是一系列的脈衝信號，那問題就轉變成了：一個信號的頻譜函數和狄拉克梳狀函數卷積，結果是什麼？——頻譜函數被周期延拓了，而且延拓的周期就是採樣頻率！

舉個例子：

對於餘弦函數而言，比如 x(t)=Acos(2picdot8cdot t) ，其傅里葉變換包含兩個頻率分量，分別是8Hz以及-8Hz，如下圖：

我們用一個20Hz的狄拉克梳狀函數（採樣函數）對餘弦信號進行抽樣，採樣函數為：

delta_s(t)=sum_{n=0}^{N}delta(t-nT_s) ,其中 T_s=1/f_s , f_s=20Hz 。

其時域信號以及頻譜為：

可見，抽樣信號的頻譜也呈現衝擊信號的樣貌，頻譜中兩個衝擊串的間隔就是採樣頻率。

採樣後的信號為：

x_s(t)=x(t)delta_s(t)=Acos(2picdot8cdot t)sum_{n=0}^{N}delta(t-nT_s)

採樣後的時域信號及頻譜為：

可見，採樣後的信號的頻譜被周期延拓了，延拓的周期就是20Hz，也就是採樣頻率。

把以上三張圖放在一起，相比能看的更清楚一些。

到這裡，對數字敏感的童鞋可能已經猜到採樣後的12Hz是哪來的了—— -8Hz平移一個採樣周期（20Hz）得來的。進而可以得到：28Hz是8Hz平移一個採樣周期來的，32Hz是-8Hz平移兩個採樣周期來的，48Hz是8Hz平移兩個採樣周期來的，....

通過以上的分析，我們得到了如下結論：對一個連續信號的採樣，採樣後的頻譜相當於將採樣前的頻譜進行了延拓，延拓的周期就是採樣頻率。

現在假設一個信號的頻譜如下：

頻譜中最大的頻率為 f_{max} ，用一個周期為 f_s 狄拉克梳狀函數進行采採樣後的頻譜為原頻譜的周期延拓，示意圖如下：

看著很完美的樣子，但是如果出現一種情況： f_{max}>frac{1}{2}f_s ，結果會怎麼樣？

也就是說，原始頻譜信號經過周期延拓後會有一部分重疊，這樣在重疊的部分就會有信息的丟失，也就無法進行信號復原了，也就是說，對於連續信號的進行抽樣離散的話，必須保證採樣頻率是原連續信號最大頻率分量的2倍頻率以上，否則信號就難以復原。這就是採樣定理，又叫奈奎斯特採樣定理或香農採樣定理。

——————————————————————————————

更多工程中的數學和物理知識，請移步潘工的專欄：

小潘工程師筆記?

zhuanlan.zhihu.com圖標
推薦閱讀：
查看原文 >>
相关文章