例說信號處理與濾波囂設計

信號分析與處理概述

博客對圖片和公式的支持不是很好，相應的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系統——出鈔機：若是你今天往出鈔機投幣口(輸入端)投幣一元，它將在接下來的三天，每一天都會從出幣口(輸出端)吐出一元。試問：若是某個月的一、三、4日，你分別投了二、二、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；濾波器結構的選擇，量化對係數的影響)