信號分析與處理概述
博客對圖片和公式的支持不是很好，相應的word版本見這里。
目錄
數字時代... 2
數字信號處理的應用... 3
頻率——信號的指紋... 5
卷積可以不卷... 8
向量運算的啟示... 11
濾波器設計征程... 16
最后一擊——濾波的實現方法... 22
縱覽全局... 27
數字時代
信號處理是對原始信號進行改變，以提取有用信息的過程，它是對信號進行變換、濾波、分析、綜合等處理過程的統稱。數字信號處理是將信號以數字方式表示並處理的理論和技術；模擬信號處理是指用模擬系統對模擬信號進行處理的方法或過程。
數字信號處理課程的主要內容包括信號分析與處理。兩者並不是孤立的，不同的信號處理方法往往需要選擇不同的信號表示形式。兩者的區別主要表現在，信號處理是用系統改變輸入信號，以得到所期望的輸出信號，如信號去噪；而信號分析往往是通過變換(傅里葉變換、小波變換等)，或其它手段提取信號的某些特征，如語音信號的基本頻率，圖像的直方圖等。
早期的信號處理局限於模擬信號，隨着數字計算機的飛速發展，信號處理的理論和方法得以飛速發展，出現了不受物理制約的純數學的加工，即算法，並確立了數字信號處理的領域。現在，對於信號的處理，人們通常是先把模擬信號變成數字信號，然后利用高效的數字信號處理器(DSP：Digital Signal Processor)或計算機對其進行數字形式的信號處理。
一般地講，數字信號處理涉及三個步驟：
(1) 模數轉換(A/D轉換)：把模擬信號變成數字信號，是一個對自變量和幅值同時進行離散化的過程，基本的理論保證是采樣定理。
(2) 數字信號處理(DSP)：包括變換域分析(如頻域變換)、數字濾波、識別、合成等。
(3) 數模轉換(D/A轉換)：把經過處理的數字信號還原為模擬信號。通常，這一步並不是必須的。
圖 1數字信號處理基本步驟
數字信號處理的應用
圖 2家庭影院(視+聽)
高分辨率圖像、高保真音質
圖 3語音識別
噪聲環境下高識別率
圖 4圖像增強
更清晰、美觀
圖 5無人駕駛
高度智能、安全
圖 6醫學成像
更高效、更精確的成像結果
狹義地講，信號處理可以統稱為濾波，根據不同的要求，選用不同性能的濾波器。在數字信號處理應用中，設計合適的濾波器至關重要。什么是數字濾波器呢？數字濾波器就是數字信號處理(Digital Signal Processing)算法，或者是數字信號處理器件(Digital Signal Processor)。什么是數字信號處理算法呢？我們需要借助基本的數學工具和方法。
首先回到信號處理的根本目的：用濾波器改變輸入信號，以期得到理想的輸出信號，如果輸入信號用$x\left[ n \right]$表示，輸出信號用$y\left[ n \right]$表示，則數字信號處理或濾波可用如下框圖表示：
圖 7數字信號濾波圖示
上述框圖並沒有給出數字信號處理算法的具體實現方法，但提供了數字信號處理的基本數學思想：通過系統將$x\left[ n \right]$和$y\left[ n \right]$進行關聯。在數學上描述$x$和$y$之間關系的手段有很多，如函數$y = f\left( x \right)$，方程組$Ax = y$等，從數學角度來理解，如果$f\left( \cdot \right)$為不同映射，即使同樣的自變量$x$，也會得到不同的因變量$y$；如果系數矩陣$A$不同，同樣的$x$經變換之后，將得到不同的$y$。在數字信號處理課程中，描述輸入輸出信號關系的基本數學工具是差分方程：
\[\sum\limits_{k = 0}^Q {{a_k}y\left[ {n - k} \right]} = \sum\limits_{k = 0}^P {{b_k}x\left[ {n - k} \right]} \]
令$N = \max \left( {P,Q} \right)$，上述差分方程可寫為
\[\sum\limits_{k = 0}^N {{a_k}y\left[ {n - k} \right]} = \sum\limits_{k = 0}^N {{b_k}x\left[ {n - k} \right]} \]
如果僅局限於基礎(經典)的數字信號處理算法研究，全部內容都是圍繞這個線性常系數差分方程(Linear Constant-Coefficient Difference Equation，LCCDE)展開的。那么，我們需要研究LCCDE什么呢？很顯然，我們可以將LCCDE看成一個系統，它描述了輸入$x\left[ n \right]$和輸出$y\left[ n \right]$之間的關系。哪些參數決定這個系統的性能呢？方程中除了$x\left[ n \right]$和$y\left[ n \right]$，還有三個參數：方程階數$N$，系數${a_k}$和${b_k}$——數字濾波器設計的根本任務就是確定這三個參數——目標非常明確!!!
從何下手呢？如何判定所設計的濾波器符合預期呢？迷茫中ing…
頻率——信號的指紋
烈日當空，窗外的知了熱得撕心裂肺地吼叫
聲音特點：又大又尖，怎叫人不心煩意亂？!
大：能量大，或者信號幅度大(聲壓級超過120dB即達到痛閾)；
尖：頻率高，即知了的叫聲中包含了豐富的高頻成分，因此聽起來很“刺耳”。
圖 8知了聲的頻譜圖(橫坐標：頻率(Hz)；縱坐標：幅度(dB))
圖 9 聽覺曲線
圖 10 聲音頻率范圍
可見，幅度和頻率是聲音信號的主要參數。諸如“低音炮”、“高音喇叭”、“男低音”和“女高音”都是從幅度和頻率去描述聲音信號的，因此，分析信號的頻率特性是信號處理領域的重要內容。
分析信號的頻率特性無非就是想知道信號包含了哪些頻率成分，各頻率成分的大小是多少。傅里葉分析完美地解決了這個問題。在信號與系統課程中，我們通常將傅里葉分析分為傅里葉級數和傅里葉變換，前者用來對付周期信號，后者用來處理非周期信號，但無論是傅里葉級數，還是傅里葉變換(當然，傅里葉級數也可以納入到傅里葉變換體系中)，它們的原理或宗旨都一樣：將一般的信號分解為基本信號的線性組合
(連續周期信號) \[\tilde x\left( t \right) = \sum\limits_{k = - \infty }^{ + \infty } {{a_k}{e^{jk{\omega _0}t}}} \]
(離散非周期序列) \[x\left[ n \right] = \frac{1}{{2\pi }}\int_{\left\langle {2\pi } \right\rangle } {X\left( {{e^{j\omega }}} \right)} {e^{j\omega n}}d\omega \]
以式為例，等式左邊是一般的(滿足收斂條件的)周期信號，而右邊則是頻率為$k{\omega _0}$、幅度為${a_k}$的虛指數信號${e^{jk{\omega _0}t}}$的線性組合。
現在可以談論一下濾波的概念了。假設式中$\tilde x\left( t \right)$就是知了發出的“嗞哇嗞哇……”沒完沒了的周期信號，現在我們想把煩人的高頻成分去掉，最直觀的想法就是需要設計一個濾波器，這個濾波器可以讓較低頻率的信號順利通過，同時又能阻止較高頻率的信號，這就是所謂的低通濾波器。
卷積可以不卷
再回到線性常系數性差分方程，參數$N$，${a_k}$和${b_k}$完全決定了方程所描述的系統的所有特性。什么？難道與輸入$x\left[ n \right]$和輸出$y\left[ n \right]$沒關系？對，沒有關系!電阻是一個系統，其阻值僅與自身的材料及結構有關，雖然有關系式$R = v\left( t \right)/i\left( t \right)$，事實上，即使兩端沒有電壓，電阻依然存在；再比如理想的線性放大電器，它的增益僅取決於內部結構，而與輸入輸出無關，但為了測量放大器的增益(放大倍數)，我們可以在輸入端接入幅度為$i$的信號，然后測量輸出信號的幅度$o$，這樣就可以得到放大倍數$g = \frac{o}{i}$，更特殊地，如果取$i = 1$，此時輸出信號的幅度就是放大器的增益，即$g = o$。類似的概念推廣到LCCDE描述的LTI系統，我們如何獲得這個系統的特性呢？輸入——輸出描述法，即令輸入$x\left[ n \right]$取某種特殊值時，計算(或測量)系統的輸出$y\left[ n \right]$。那么，什么樣的$x\left[ n \right]$算特殊呢？單位脈沖!
\[\delta \left[ n \right] = \left\{ \begin{array}{l}
1,\;\;\;n = 0\\
0,\;\;n \ne 0
\end{array} \right.\]
也就是說，令$x\left[ n \right] = \delta \left[ n \right]$，此時系統的輸出就是所謂的單位脈沖響應$h\left[ n \right]$。再與放大器的例子對比一下：輸入信號幅度$i = 1$，輸出信號的幅度$o$就是放大倍數(放大器的重要指標，從應用者角度而言，其實就是唯一關心的指標)。我們有理由相信，對於LCCDE描述的LTI系統，得到了單位脈沖響應$h\left[ n \right]$，就能夠掌握系統的全部特性(因果性、穩定性、頻率選擇性等)。何以見得呢？之后將逐步解答。
現在，我們假設已經知道了LTI系統的單位脈沖響應$h\left[ n \right]$，對於任意輸入$x\left[ n \right]$，如何求系統的輸出$y\left[ n \right]$呢？其實就是求濾波后的信號。首先，我們要建立已知和需求之間的聯系。已知的是 ，需求是 。
圖 11 單位脈沖響應
事實上，任意序列$x\left[ n \right]$很容易用$\delta \left[ n \right]$的移位、加權的線性組合來表示：$x\left[ n \right] = \sum\limits_{k = - \infty }^{ + \infty } {x\left[ k \right]\delta \left[ {n - k} \right]} $。
圖 12 序列的單位脈沖表示
結合系統的線性和時不變性，有
式中最后一個等式\[y\left[ n \right] = \sum\limits_{k = - \infty }^{ + \infty } {x\left[ k \right]h\left[ {n - k} \right]} \]就是卷積和，記作$y\left[ n \right] = x\left[ n \right]*h\left[ n \right]$，表明LTI系統的(零狀態)響應$y\left[ n \right]$是單位脈沖響應$h\left[ n \right]$的移位、加權、求和。這句話也描述了求$y\left[ n \right]$的步驟。根據圖 11，我們可以假想有這樣的一個LTI系統——出鈔機：如果你今天往出鈔機投幣口(輸入端)投幣一元，它將在接下來的三天，每一天都會從出幣口(輸出端)吐出一元。試問：如果某個月的1、3、4日，你分別投了2、2、3元，那么，出鈔機每天的輸出是多少元呢？不妨用圖來說明一下。
圖 13 卷積和的圖解
圖解結果告訴我們，出鈔機在2-7日分別吐出2,2,4,5,5,3元。圖解的過程包含了乘積與求和(“積”與“和”)，但並有體現出“卷”(翻轉)這一操作。如果只想求某一天(比如5日)出鈔口吐出多少錢，此時就要用另一種方法，即許多教材中描述的步驟：
(1) 變量替換：$x\left[ n \right] \to x\left[ k \right],h\left[ n \right] \to h\left[ k \right]$
(2) 將$x\left[ k \right]$或$y\left[ k \right]$翻轉
(3) 移位、相乘、相加
以上步驟包含了“卷積和”所有的操作(卷——翻轉；積——相乘；和——相加)。
圖解過程是根據線性時不變系統的定義導出的一種結果，是系統特性的直接反映。卷積和是信號與系統、數字信號處理中最重要的公式之一，它描述了LTI系統輸入$x\left[ n \right]$、輸出$y\left[ n \right]$以及單位脈沖響應$h\left[ n \right]$之間的關系，已知三者中的任何兩個，就可以確定第三者，於是就有以下應用場合：
(1) 已知輸入$x\left[ n \right]$和系統單位脈沖響應$h\left[ n \right]$，求輸出$y\left[ n \right]$，即用確定的系統對輸入信號濾波處理；
(2) 已知輸出$y\left[ n \right]$和系統單位脈沖響應$h\left[ n \right]$，求輸入$x\left[ n \right]$；
(3) 已知輸入$x\left[ n \right]$和輸出$y\left[ n \right]$，求系統單位脈沖響應$h\left[ n \right]$。
現在，我們將注意力集中到卷積和公式上來，公式重寫如下
\[y\left[ n \right] = \sum\limits_{k = - \infty }^{ + \infty } {x\left[ k \right]h\left[ {n - k} \right]} \]
卷積和只包含了移位、數乘和相加運算，看來，數字信號處理的計算十分簡單。如果再上升一個層次，卷積和能很直觀解釋系統如何改變輸入信號的嗎？從圖 13的結果很難看出系統是如何將$x\left[ n \right]$變成$y\left[ n \right]$的，而且也很難說明這種變化意味着什么，因此，卷積和在解釋系統行為特性時很不直觀，這是卷積和運算本身決定的。
向量運算的啟示
了解運算背后的物理或幾何含義對理解問題的本質至關重要。例如向量的乘法和加法：設兩個向量$\vec A = {x_1} + i{y_1}$，$\vec B = {x_2} + i{y_2}$，向量的加法符合平行四邊形法則，從圖 14很容易看出兩個向量相加的幾何意義，而相應的代數運算也很簡單——對應分量相加即可，即有
$\vec C = \vec A + \vec B = \left( {{x_1} + {x_2}} \right) + i\left( {{y_1} + {y_2}} \right)$
但如果要求兩個向量的乘法呢？代數運算的結果是
相比兩個向量的加法，向量的乘法運算要復雜一些，更重要的是，根據運算結果，我們很難看清向量乘法的幾何意義。如果將向量$\vec A$和$\vec B$用極坐標表示，情況又會如何呢？設
則兩個向量的乘法
計算簡單，幾何意義明確——模相乘，相位相加。
圖 14 向量加的直角坐標表示
圖 15 向量乘法的極坐標表示
向量運算的例子告訴我們，選擇不同的坐標系，不僅影響運算的復雜度，在解釋運算的幾何意義時也各有千秋。上例告訴我們，如果是向量的加法，直角坐標系比極坐標系方便；如果是向量的乘法，很顯然，極坐標系比直角坐標系方便。
回到數字信號處理的話題，卷積和是讓無數人困惑的公式，然而，它又是經典數字信號處理算法的根源。既然卷積和運算就是濾波，我們如何評判濾波的效果呢？從什么角度理解系統的濾波行為更好呢？
再想想知了的叫聲，之所以煩人，是因為包含了大量的高頻分量。看來，用頻率這一參數來解釋信號特性很符合人的直觀感覺。我們有更好的方式來洞悉卷積和的物理意義嗎？能不能站在頻(率)域的角度來解釋卷積和呢？當然可以，因為我們有卷積定理(性質)：
通過傅里葉變換，將時域中的卷積和運算換成了頻域中的相乘運算。那么頻域中的相乘運算有什么好處呢？不要忘了，我們主要是為了直觀解釋濾波的行為和特性。
圖 16是一段濾波前知了聲$x\left[ n \right]$的時域波形，單從圖形上看，似乎很難說出這一段信號的特點，運用傅里葉變換，得到$x\left[ n \right]$的幅度譜$\left| {X\left( {{e^{j\omega }}} \right)} \right|$如圖 17所示，很明顯，該圖說明信號$x\left[ n \right]$包含了豐富的高頻成分。現在，我們希望設計一個濾波器，可以濾除2500Hz以上的頻率分量，目的是僅保留紅框內的頻率分量，因此要求濾波器的頻率響應(幅頻響應)是
濾波器的幅頻特性如圖 18所示，截止頻率${f_c} = 2500$Hz。由卷積定理可知，輸出信號的幅度譜為
\[\left| {Y\left( {{e^{j\omega }}} \right)} \right| = \left| {X\left( {{e^{j\omega }}} \right)} \right|\left| {H\left( {{e^{j\omega }}} \right)} \right|\]
從圖形上來看就更直觀了：圖 17與圖 18相乘，得到了濾波后信號的幅度譜，濾波前、后信號頻譜的變化一目了然，很顯然，它是按照濾波器頻率響應$H\left( {{e^{j\omega }}} \right)$所規定的形狀去改變的。而對比濾波前、后的時域波形，很難解釋信號的什么特性被改變了。
在上述的討論中，我們只關注幅度譜，實際中，相位譜也必須考慮，更完整的表達式為
對比極坐標系下兩個向量的乘法，能找到相似之處嗎？——模相乘、相位相加。如果將卷積和$y\left[ n \right] = x\left[ n \right]*h\left[ n \right]$看作直角坐標系下兩個向量的乘法，那么\[Y\left( {{e^{j\omega }}} \right) = X\left( {{e^{j\omega }}} \right)H\left( {{e^{j\omega }}} \right)\]就類似於極坐標系下兩個向量的乘法——運算簡單、物理意義直觀。[通過向量乘法運算的例子，我們可以看到，向量的表示方式起到了關鍵性的作用；同樣地，信號也有不同的描述方式，如時域、頻域、復頻域，不同的表示方法都是為了使分析問題更簡單、物理意義更明確]
圖 16濾波前知了聲$x\left[ n \right]$時域波形
圖 17濾波前知了聲的幅度譜$\left| {X\left( {{e^{j\omega }}} \right)} \right|$ (紅框是希望保留的低頻分量)
圖 18低通濾波器幅頻特性$\left| {H\left( {{e^{j\omega }}} \right)} \right|$ (藍色線)
圖 19濾波后知了聲的幅度譜$\left| {Y\left( {{e^{j\omega }}} \right)} \right|$
圖 20濾波后知了聲$y\left[ n \right]$時域波形
濾波器設計征程
借助系統的頻率響應$H\left( {{e^{j\omega }}} \right)$，可以很方便地解釋濾波器的特性(主要是指頻率選擇性)。有什么定量指標來描述頻率選擇性嗎？通常有四個指標——橫坐標兩個，縱坐標兩個。以低通濾波器為例來說明。如圖 21所示，橫坐標兩個指標為通帶截止頻率${\omega _p}$和阻帶截止頻率${\omega _s}$，縱坐標兩個指標為通帶衰減${A_p}$和阻帶衰減${A_s}$。
圖 21數字低通濾波器指標
通常而言，濾波器設計就是指根據應用的需求，提出合理的濾波器指標(通帶截止頻率和衰減，阻帶截止頻率和衰減)，然后借助濾波器設計工具得到濾波器參數(系數)。濾波器系數？沒聽說過啊，是什么東西？其實我們早見過面了，就是差分方程中的${a_k}$和${b_k}$!早就說過，數字濾波器設計要不忘初心，根本任務就是確定這三個參數——方程的階數$N$，系數${a_k}$和${b_k}$。
如何建立起濾波器指標和濾波器系數之間的聯系呢？很顯然，濾波器的指標約束着其頻率響應$H\left( {{e^{j\omega }}} \right)$，而$H\left( {{e^{j\omega }}} \right)$與LCCDE的階數$N$，系數${a_k}$和${b_k}$直接相關。傅里葉變換又該出場了。對差分方程兩邊取傅里葉變換得
\[\sum\limits_{k = 0}^N {{a_k}{e^{ - jk\omega }}Y\left( {{e^{j\omega }}} \right)} = \sum\limits_{k = 0}^N {{b_k}{e^{ - jk\omega }}X\left( {{e^{j\omega }}} \right)} \]
從而可得頻率響應
\[H\left( {{e^{j\omega }}} \right) = \frac{{Y\left( {{e^{j\omega }}} \right)}}{{X\left( {{e^{j\omega }}} \right)}} = \frac{{\sum\limits_{k = 0}^N {{b_k}{e^{ - jk\omega }}} }}{{\sum\limits_{k = 0}^N {{a_k}{e^{ - jk\omega }}} }}\]
易見，頻率響應$H\left( {{e^{j\omega }}} \right)$完全由$N$，${a_k}$和${b_k}$確定。綜上可知，濾波器設計就是如何根據濾波器指標求解濾波器系數的過程，從數學角度來看，就是函數逼近和數值優化的過程。經典濾波器設計理論發展較完備(巴特沃斯、切比雪夫I、切比雪夫II、橢圓……，各自的特點是什么？更多內容，請點這里)，許多教科書都有詳細討論，此處不再贅述，我們的重點是要借助現代化的設計工具，輕松完成設計任務並實施濾波處理。此處僅以MATLAB工具為例，並通過兩種方式來實現濾波器設計——腳本代碼和可視化方法。
需求：設計一個低通濾波器來處理知了聲，通帶截止頻率${F_p} = 2000$Hz，通帶衰減${A_p} = 0.1$dB，阻帶截止頻率${F_s} = 2500$Hz，通帶衰減${A_p} = 50$dB。
乍一看，這些指標與圖 21中的指標完全不同，實際需求的指標往往同本例，頻率指標為模擬頻率，衰減指標以dB為單位，而設計數字濾波器時，往往需要數字頻率指標。如何進行指標轉換呢？直接用公式計算${\omega _p} = 2\pi {F_p} = 4000\pi $，${\omega _s} = 2\pi {F_s} = 5000\pi $嗎？圖 21告訴我們，無論通帶截止頻率，還是阻帶截止頻率都不會超過$\pi $。難道最大的數字角頻率就是$\pi $嗎？為什么？這個問題同樣令許多人困惑。如何衡量數字頻率的大小呢？先考慮一個簡單的數字信號：$x\left[ n \right] = \cos \left( {\omega n} \right)$，針對$\omega $取不同的值，分別畫出對應的波形如下：
圖 22$\cos \left( {\omega n} \right)$波形圖
振盪得越厲害，意味着頻率越高，這是符合常理的。上圖中的9個波形，哪個振盪得最厲害呢？第二排中間那個——+1和-1交替，而此時$\omega = \pi $。而且還發現，$\omega = 0$和$\omega = 2\pi $的波形一樣，並不是$\omega $的值越大，振盪得越快。如果是連續時間信號$x\left( t \right) = \cos \left( {2\pi ft} \right) = \cos \left( {\Omega t} \right)$，大家都知道，隨着$\Omega $的增大，波形會振盪得越快。從連續時間信號到離散時間信號，是什么原因導致兩者出現如此大的差異呢？通過對采樣過程的探索，我們便可理解其中的奧妙。設以時間間隔${T_s}$對$x\left( t \right) = \cos \left( {\Omega t} \right)$采樣，這樣便得到
$x\left( {n{T_s}} \right) = \cos \left( {\Omega n{T_s}} \right) = \cos \left( {\omega n} \right)$
因此有
$\omega = \Omega {T_s} = \Omega /{f_s}$
其中${f_s} = 1/{T_s}$為采樣頻率，$\Omega $為模擬角頻率，$\omega $為數字角頻率。結論是：數字頻率是模擬頻率關於采樣頻率的歸一化。
根據Nyquist采樣定理，當采樣頻率大於信號最高頻率2倍時，可以由樣點完全恢復原來的信號。我們不妨假設信號的最高頻率正好是采樣頻率的一半，即\[{f_{\max }} = {f_s}/2\]，對應的最高模擬角頻率${\Omega _{\max }} = 2\pi {f_{\max }} = \pi {f_s}$，代入式得
${\omega _{\max }} = {\Omega _{\max }}/{f_s} = \pi {f_s}/{f_s} = \pi $
這也證明了為什么最高數字角頻率為$\pi $。(如果不滿足Nyquist采樣條件又會怎樣呢？事實上並不影響結論成立)。關於各種頻率之間的關系的討論，可以參考這里(注意文中“fdtool工具中歸一化頻率為什么最大只到0.5的原因”，應改為“最大只到1的原因”。)
再回到濾波器設計題目，第一步，需要將模擬頻率指標轉換為歸一化數字頻率指標，轉換公式是
\[\omega = 2f/{f_s}\;\;\;\;(\pi )\]
很顯然，我們需要知道采樣頻率${f_s}$。知了聲的采樣頻率${f_s} = 11025$Hz。將模擬頻率指標代入式計算出對應的歸一化數字頻率指標
這樣，我們就得到了數字域的四個指標：${\omega _p},{A_p},{\omega _s},{A_s}$。如何從這四個指標去計算濾波器系數呢？有現成的公式套用，即許多教科書上所說的各類濾波器原型。我們選擇橢圓型濾波器為例。MATLAB代碼如下：
clc;
close all
fs=11025; %采樣頻率(Hz)
Fp = 2000; %通帶截止頻率(Hz)
Fs = 2500; %阻帶截止頻率(Hz)
Ap = 0.1; %通帶衰減(dB)
As = 50; %阻帶衰減(dB)
wp=2*Fp/fs; %歸一化通帶截止頻率
ws=2*Fs/fs; %歸一化阻帶截止頻率
[N,Wp]=ellipord(wp,ws,Ap,As); %確定帶通濾波器的階數和截止頻率
[b,a]=ellip(N,Ap,As,Wp); %確定濾波器系數
[h,w]=freqz(b,a); %求數字帶通濾波器的頻率響應
%以下為繪圖命令，繪制帶通濾波器的幅頻響應
figure;
plot(w*fs/(2*pi),20*log10(abs(h)/max(abs(h))));
axis([0,fs/2,-100,0]);
title('數字低通濾波器的幅度響應');
xlabel('頻率(Hz)');
ylabel('幅度(dB)');
grid
濾波器系數
a = [1 -2.98876196757509 5.41154092244290 -6.18465137960809 4.94650037623018 -2.65349822687766 0.903727225595414 -0.150115115698030]
b = [0.0204237714181223 0.0185281305134523 0.0518202070683316 0.0515988082549051 0.0515988082549052 0.0518202070683315 0.0185281305134523 0.0204237714181222]
分別對應式中的${a_k}$和${b_k}$，得到了這兩個參數，濾波器設計基本上算完成了。
如果你掌握了濾波器設計的基本理論，但又不想寫代碼，你可以利用FDATool輕松完成濾波器設計的全過程。
在命令窗口輸入：fdatool，回車，就可以看到下圖所示的界面：
然后動動鼠標就可以完成濾波器設計了。針對本例，我們按下圖紅框所示設定指標，在完成指標設定之后，點“Design Filter”按鈕即完成設計。關於濾波器的所有信息，都可以通過菜單或工具欄獲得。注意”Current Filter Information”部分，它描述了當前所設計出的濾波器基本信息，包括結構、階次、穩定性等。有一個問題值得思考：濾波器的結構對濾波器的工程實現有什么影響呢？成本、穩定性、抗噪聲性能，這些都與結構有關，因此，掌握濾波器各種結構的優缺點也十分重要。
很顯然，用FDATool設計濾波器十分簡單!無論采用什么方法，都不要忘記濾波器設計的根本目標——得到濾波器系數。關於FDATool更多內容，請點這里和這里。
當濾波器設計完成之后，我們該考慮如何用所得到的系數對輸入信號進行濾波了。
最后一擊——濾波的實現方法
當濾波器系數${a_k}$和${b_k}$確定之后，方程也就完全確定了，所謂的濾波，就是要對給定的輸入$x\left[ n \right]$，求系統的響應$y\left[ n \right]$。MATLAB提供了幾種濾波實現函數：
conv：線性卷積，可以用來實現濾波
filter：濾波，但是輸出點數和輸入點數相同，其余的存在內部狀態中，以便處理后續輸入
filtfilt：零相位濾波，輸出通過濾波器后再反向通過濾波器，消除相位影響
fftfilt：利用基於FFT的重疊相加法對數據進行濾波，這種頻域濾波技術只對FIR濾波器有效
完整例子：用MATLAB讀入知了聲的音頻文件，畫時域波形和頻譜，設計濾波器，對音頻文件進行濾波，最后畫濾波后的信號時域波形和頻譜。程序如下：
clc;
close all
[xn, fs] = wavread('cicada.wav'); %讀wav文件(知了聲)，獲得樣點值xn和采樣頻率fs
%計算信號xn的頻譜
xn_FFT = fft(xn); %FFT頻譜分析
xn_FFT_AMP = abs(xn_FFT); %求xn的幅度譜
xn_FFT_AMP_db = 20*log10(xn_FFT_AMP/max(xn_FFT_AMP)); %用對數表示
%畫xn時域波形
figure
n = 1:length(xn); %生成樣點橫坐標向量
plot(n,xn); %以n為橫坐標，畫xn
axis([1 length(xn) -1 1]) %設定坐標軸顯示范圍
xlabel('樣點'); %設定x坐標標記
ylabel('幅度'); %設定y坐標標記
title('知了聲時域波形'); %設定標題
%畫xn的幅度譜(對數形式)
figure
faxis = (0:length(xn)-1)*fs/length(xn); %生成橫坐標(頻率軸)刻度，FFT的頻率范圍[0 fs]
plot(faxis,xn_FFT_AMP_db); %以faxis為橫軸畫xn頻譜圖
axis([0 fs/2 -100 0]) %設定顯示范圍，由於FFT頻譜具有對稱性，因此只顯示前半部分
xlabel('頻率(Hz)');
ylabel('幅度(dB)');
title('知了聲幅度譜');
% 設計低通濾波器
Fp = 2000; %通帶截止頻率(Hz)
Fs = 2500; %阻帶截止頻率(Hz)
Ap = 0.1; %通帶衰減(dB)
As = 50; %阻帶衰減(dB)
wp=2*Fp/fs; %歸一化通帶截止頻率
ws=2*Fs/fs; %歸一化阻帶截止頻率
[N,Wp]=ellipord(wp,ws,Ap,As); %確定帶通濾波器的階數和邊緣頻率
[b,a]=ellip(N,Ap,As,Wp); %確定濾波器系數
[h,w]=freqz(b,a); %求數字帶通濾波器的頻率響應
%以下為繪圖命令，繪制帶通濾波器的幅頻響應
figure;
plot(w*fs/(2*pi),20*log10(abs(h)/max(abs(h))));
axis([0,fs/2,-100,0]);
title('數字低通濾波器的幅度響應');
xlabel('頻率(Hz)');
ylabel('幅度(dB)');
grid
%利用設計的濾波器對xn進行濾波，得到濾波輸出yn
yn = filter(b,a,xn); %濾波函數，b和a為濾波器系數，xn是輸入信號
wavwrite(yn,fs,8,'cicada_LPFiltered.wav'); %保存濾波后的信號
%畫濾波后的信號yn時域波形
figure
n = 1:length(yn);
plot(n,yn);
axis([1 length(yn) -1 1])
xlabel('樣點');
ylabel('幅度');
title('低通濾波后知了聲時域波形');
%計算濾波后信號yn的頻譜
yn_FFT = fft(yn);
yn_FFT_AMP = abs(yn_FFT);
yn_FFT_AMP_db = 20*log10(yn_FFT_AMP/max(yn_FFT_AMP));
figure
faxis = (0:length(yn)-1)*fs/length(yn);
plot(faxis,yn_FFT_AMP_db);
axis([0 fs/2 -100 0])
xlabel('頻率(Hz)');
ylabel('幅度(dB)');
title('濾波后知了聲幅度譜');
結果：
圖 23知了聲時域波形
圖 24知了聲頻譜(幅度頻)
圖 25低通濾波器頻率響應(幅頻響應)
圖 26低通濾波后知了聲時域波形
圖 27濾波后知了聲頻譜圖(幅度譜)
縱覽全局
經典信號處理課程的主要內容：
(1) 信號的表示方法(時域表示、頻域表示、復頻域表示)
(2) 系統的描述方法(LCCDE、單位脈沖響應、頻率響應、系統函數、零極點圖)
(3) 濾波器設計方法(代碼和可視化方法)
(4) 濾波的實施方法(濾波器類型的選擇——IIR、FIR；濾波器結構的選擇，量化對系數的影響)