《數(shù)字信號處理教程-MATLAB釋義與實現(xiàn)》第九章

第九章

綜合應(yīng)用1綜合應(yīng)用本章內(nèi)容要點：9.1信號處理工具（sptool）的介紹9.2語音數(shù)字化量化噪聲的改善9.3系數(shù)量化和運算量化的影響9.4 在雙音撥號系統(tǒng)中的應(yīng)用9.5正余弦信號的譜分析

9.6音樂信號處理9.7變采樣率數(shù)字濾波9.8稀疏天線陣列設(shè)計9.9結(jié)束語2綜合應(yīng)用信號處理的概念比較深，其中許多概念不是單純用數(shù)學(xué)推導(dǎo)就能掌握的，往往要通過在實踐中反覆運用，才能真正掌握。開辟這一章的目的就是通過一些實例來說明如何把理論用到工程實踐中去。同時又加深了對理論的理解。另外在本章中還要介紹一些MATLAB的工具，本章先介紹信號處理工具sptool，它是信號處理工具箱中的一個集成環(huán)境。學(xué)習(xí)和使用這個工具可以幫助讀者把已學(xué)的信號處理知識系統(tǒng)化和集成化，對于在工程中進(jìn)行信號處理會有很大的幫助。3綜合應(yīng)用

9.2節(jié)介紹信號的量化問題，將介紹用壓縮擴(kuò)張器減小相對誤差的非線性處理方法；9.3節(jié)介紹濾波器系數(shù)量化和計算量化的影響；9.4節(jié)介紹雙頻撥號系統(tǒng)，它是一個從雙頻發(fā)送到接收檢測的完整的系統(tǒng)，其中包括goertzel算法開發(fā)；9.5節(jié)介紹了頻譜分析儀中的數(shù)字信號處理技術(shù)；9.6節(jié)介紹音響系統(tǒng)中的信號處理，它涉及回聲和混響的生成和處理；9.7節(jié)介紹變采樣率系統(tǒng)，包括內(nèi)插和抽取，這也屬于一種非線性處理；9.8節(jié)把時域信號擴(kuò)展到空域信號處理；4綜合應(yīng)用上述的七節(jié)都是以前面學(xué)過的八章為基礎(chǔ)，但都把原有的知識拓寬和加深了一步。在本書中，我們對每一個問題都只作了簡要的提示，各個問題都有很大的展開余地，所以其中任何一個問題都可以作為實驗或課程設(shè)計的基礎(chǔ)。這要由教師和讀者自行取舍，進(jìn)行適當(dāng)?shù)拈_發(fā)。9.9節(jié)是結(jié)束語，大略地介紹了一下數(shù)字信號處理器的應(yīng)用領(lǐng)域，并介紹了MATLAB中與信號處理有關(guān)的各種工具箱函數(shù)。59.1信號處理工具介紹

(講本節(jié)時，應(yīng)盡量直接用MATLAB的圖形界面)信號處理工具sptool（SignalProcessingTool的縮寫）為信號處理的研究工作提供了一個集成環(huán)境和工具。信號處理的任務(wù)：一是對信號進(jìn)行分析；二是濾波器設(shè)計。在這兩個任務(wù)進(jìn)行過程中，經(jīng)常要做第三個任務(wù)，那就是要把信號加到所設(shè)計的濾波器中，看它的輸出是否滿足要求。Sptool把這三個任務(wù)集成在一起，加上了適當(dāng)?shù)墓芾砉δ?，配以良好的工作界面，利用它可以大大提高研究工作的效率?/p>

6信號處理工具(sptool)介紹 在MATLAB命令窗中，鍵入fdatool,得到如右圖的界面。各欄中分別存入了系統(tǒng)中原來已保存的信號、濾波器和頻譜的名稱。7信號處理工具(sptool)介紹1．

信號和濾波器的導(dǎo)入先在MATLAB工作空間中放入所需的分析的信號和濾波器。設(shè)濾波器就是例8.8.1中設(shè)計并導(dǎo)入了工作空間的濾波器，再來建立一組信號。在命令窗中，鍵入：n=0:200;T=1;s1=sin(5*2*pi*n*T);s2=sin(10*2*pi*n*T);s3=cos(15*2*pi*n*T);s=s1+s2+s3;這時s1,s2,s3和它們的合成信號s都已輸入工作空間。8信號處理工具(sptool)介紹點擊左上角的【File】及下拉菜單的【Import】,出現(xiàn)如下圖的新視窗。它分成三欄；左邊一欄【Source】，指定數(shù)據(jù)來源；如果來源是工作空間，中間欄的【W(wǎng)orkspaceContents】將顯示工作空間中的全部變量；右邊一欄指定導(dǎo)入的目標(biāo)，即要說明導(dǎo)入的數(shù)據(jù)是作為信號、還是濾波器還是頻譜。右欄的下部有【SamplingFrequency】和【Name】兩個框需要用戶填寫。按照說明的方法，可依次將信號s和8.6節(jié)設(shè)計的濾波器輸入sptool。

9信號處理工具(sptool)介紹10信號處理工具(sptool)介紹2.信號的時域和頻域觀測在圖9.1.1上，選定sig7,點擊【view】(觀測),就進(jìn)入【SignalBrowser】(信號瀏覽器)視窗。如圖9.1.3所示。要想求得該信號的頻譜，也先選定sig7，然后在頻譜欄下點擊【creat】,此時出現(xiàn)頻譜觀測器的界面，見圖9.1.4。先要在左欄上部選定求頻譜的方法，目前我們只學(xué)了一種，即選FFT.在下面的框中填入點數(shù)，例如1024，然后點擊【Apply】,就出現(xiàn)了頻譜的曲線。11信號處理工具(sptool)介紹12信號處理工具(sptool)介紹13信號處理工具(sptool)介紹3．

濾波器的觀測和修改選定濾波器的名稱filt7,點擊其下方的【View】，就可以觀測它的頻率響應(yīng)。如點擊【Edit】，那就不僅能觀測，還能修改。4.

讓信號通過濾波器求輸出信號在圖9.1.1中，選定信號和濾波器名稱，點擊濾波器欄下的按鈕【Apply】,這意味著將所選信號加到所選濾波器中去。此時將產(chǎn)生一個小視窗，提示用戶輸入信號和濾波器的名稱，并提供輸出信號缺省名稱sig8。如無修改，點擊OK，輸出信號將自動生成并列寫在信號欄中。14信號處理工具(sptool)介紹在本例中濾波器輸出sig8基本上是一個較純的10Hz正弦波，這是不難想象的。經(jīng)過仔細(xì)觀測和分析，如不滿意，可以修改濾波器后重新實驗。用這個工具，不難在很短的時間內(nèi)完成多個濾波器的設(shè)計，并分別觀測它們的輸出。所有的結(jié)果都存儲在案，可用于試驗報告。所以這是一個提高效率的有效工具。159.2語音量化噪聲的改善脈沖編碼調(diào)制(PulseCodeModulation—PCM)是把模擬信號量化為二進(jìn)制數(shù)的最簡單的方法。以N個脈沖表示N位二進(jìn)制數(shù)，以脈沖的有無判斷它是0或1。它也是用數(shù)字方式傳輸或存貯信號的常用方法之一。PCM被廣泛應(yīng)用于電話通信和利用無線電傳輸?shù)倪b測系統(tǒng)中。通過電話線傳輸?shù)恼Z音信號頻帶限于4kHz以下。因此其采樣頻率取8KHz(樣本數(shù)/秒)，并用N位二進(jìn)制序列表示它的值，每個樣本量化為2N個電平之一。所以，傳輸數(shù)字化語音信號所要求的速率為每秒8000×N位。16語音量化噪聲的改善量化處理的數(shù)學(xué)模型為

xq(n)=x(n)+q(n) (9.2.1)其中xq(n)表示x(n)的量化值，q(n)表示量化誤差，將其看作一加型噪聲。假設(shè)采用的是均勻的量化器，則可用如下均勻概率密度函數(shù)p(q)統(tǒng)計描述量化噪聲特性: (9.2.2)其中，量化步長為△=2–N。量化誤差的均方值為: (9.2.3) 17語音量化噪聲的改善用分貝來度量的噪聲均方值為: (9.2.4)可以看出，上述的量化器每增加一位，量化噪聲減小6dB，高質(zhì)量語音要求每個樣本至少量化為12位，因此傳送速率至少為96000位/秒。最大幅度為±V伏的N位（不含符號位）二進(jìn)制A/D變換器的數(shù)學(xué)模型建立如下。它把電壓V分解為2N-1份，故量化步長為V/(2N-1)，得出二進(jìn)制量化子程序bqtize。

18語音量化的子程序functiony=bqtize(x,N,V)ifnargin<3V=max(abs(x));end %V缺省取最大xax=abs(x);%去掉符號deltax=V/(2^N-1);%求量化步長xint=fix(ax./deltax+0.5); %x的量化整數(shù)y=sign(x).*xint.*deltax;%恢復(fù)量化原值這個A/D變換子程序的輸入是連續(xù)模擬電壓x，輸出則是量化了的模擬電壓y。均勻量化器在信號的整個動態(tài)范圍中的量化步長相同，所以量化噪聲均方值不變。

19語音量化噪聲的改善均勻量化器在信號的整個動態(tài)范圍中的量化步長相同，所以量化噪聲均方值不變。然而，語音信號的特性是小幅度比大幅度出現(xiàn)得頻繁。對小信號而言，量化噪聲使信噪比大大下降。解決的途徑之一是用非均勻量化器。不過在技術(shù)上制造非均勻量化器的芯片是困難的。得到非均勻量化器特性的另一個方法是用壓縮—擴(kuò)張器。可先使信號通過壓縮幅度的非線性器件，后面再接一均勻量化器，再用逆向擴(kuò)張幅度的非線性器件恢復(fù)信號。如下圖。

20語音量化的壓縮擴(kuò)張器壓縮器的作用是把小信號放大，大信號縮小，所以把壓縮器后的信號進(jìn)行量化，小信號的信噪比就得到提高。而后通過擴(kuò)張器，信號和量化噪聲同時作非線性變換，信號復(fù)原為原來電平，小信號信噪比則保持較小的水平。21語音量化的壓縮擴(kuò)張器在我國的通信系統(tǒng)中使用的對數(shù)壓縮器(稱為μ律壓縮器)具有如下輸入輸出幅度特性:其中，x是歸一化輸入，y是歸一化輸出，sign(.)是符號函數(shù)，μ=255,是控制壓縮特性的參數(shù)。由量化數(shù)字信號恢復(fù)模擬信號時，逆關(guān)系μ律為: 22語音量化的壓縮擴(kuò)張器可以用MATLAB把壓縮器和擴(kuò)張器用函數(shù)程序mulawcom和mulawexp來表示，例如把壓縮器函數(shù)程序列寫如下：functiony=mulawcom(x,mu,V)ifnargin<3|V<max(abs(x))V=max(abs(x));endy=V/log(1+mu)*log(1+mu/V*abs(x)).*sign(x);對μ律擴(kuò)張器和A-律的壓擴(kuò)器也可以按類似的方法寫出相應(yīng)的子程序。23語音量化壓縮器的特性例9.2.1畫出μ律和A律的壓縮器輸入輸出曲線。解：寫出如下的MATLAB程序hc921x=0:0.01:1;y=mulawcom(x,255,1);y1=Alawcom(x,87.56,1);plot(x,y,x,y1,':');gridonlegend('\mu律','A律')運行此程序的結(jié)果見圖9.2.1。請注意二者非常相似。24語音量化壓縮器的特性25語音量化噪聲的改善例9.2.2 設(shè)某A/D變換器把最大輸入為5伏的信號量化為四位二進(jìn)制（不含符號位），要求用圖形描寫其輸入輸出關(guān)系，并畫出其絕對誤差和相對誤差的曲線。又若信號像圖9.2.1那樣經(jīng)過壓縮擴(kuò)張器，則輸入輸出關(guān)系有何變化？解：A/D變換器的數(shù)學(xué)模型為絕對量化函數(shù)bqtize，已經(jīng)在前面得出。利用這個函數(shù)，加上求絕對誤差和相對誤差的語句，可以方便地列出以下的MATLAB程序hc922。26語音量化的信噪特性對比x=-5:0.01:5;%輸入自變量數(shù)組xq=bqtize(x,4,5);%求量化輸出e=x-xq;er=e./abs(x);%求相對誤差erplotyy(x,xq,x,er)%畫輸出及誤差曲線plot(x,e,'-.','linewidth',3)運行此程序所得曲線見圖9.2.3左圖。如在第二條語句的前后分別加上壓擴(kuò)語句：x1=mulawcom(x,255,V);%信號經(jīng)過壓縮器yq=mulawexp(xq,255,V);%信號經(jīng)過擴(kuò)張器則所得曲線見圖9.2.3右圖。27語音量化的信噪特性對比28語音量化的信噪特性對比從左圖中可以看出，在輸入為-5~5V范圍內(nèi)，絕對量化誤差e呈等幅鋸齒波形式，其最大值恒定，因而其相對誤差er在小的輸入幅度x處急遽增大，超過了1。

由右圖從中可以看出它的絕對誤差隨x的增加而增大，其相對誤差則在整個輸入范圍內(nèi)呈等幅鋸齒波形式，最大不到0.2。因此量化誤差yq-x1引起的信噪比將比不用壓擴(kuò)器有顯著的提高。299.3系數(shù)和運算量化的影響要用計算機(jī)處理信號，就不可避免地要遇到量化問題。這里包括三方面:(1)信號的量化；(2)系統(tǒng)參數(shù)的量化；以及(3)信號與系統(tǒng)相互運算中的量化。在MATLAB中，參數(shù)和運算都有十六位十進(jìn)制（64位二進(jìn)制浮點數(shù)）的精度，量化誤差通常可以忽略不計。當(dāng)要把用MATLAB設(shè)計的濾波器付諸實現(xiàn)，特別是用于實時系統(tǒng)時，要采用嵌入式的數(shù)字信號處理芯片。它們通常只有較短的（8位、16位、32位）二進(jìn)制定點數(shù)精度，必須考慮量化的影響。絕不可盲目輕信MATLAB中計算出的結(jié)果。30系數(shù)和運算量化的影響MATLAB按二進(jìn)制雙精度格式來表示數(shù)，表達(dá)一個數(shù)需要8個字節(jié)，也就是64個二進(jìn)制位。按照IEEE標(biāo)準(zhǔn)，雙精度浮點數(shù)η用下式表示其中M是一個小于一大于1/2的二進(jìn)制分?jǐn)?shù)，稱為尾數(shù)，占用52個二進(jìn)制位，近似16位十進(jìn)制；而指數(shù)E是一個帶符號的二進(jìn)制整數(shù)。占11個二進(jìn)制位，總共可表示2048個整數(shù)，即可以表示從-1023到+1024的數(shù)集，它決定了數(shù)的動態(tài)范圍。數(shù)的正負(fù)號反映在S上，它只占一位。所以用MATLAB計算通常不必考慮量化誤差。31系數(shù)和運算量化的影響系數(shù)量化：系數(shù)和信號的不同在于它的值并非來自A/D變換，而是由設(shè)計者在計算機(jī)內(nèi)設(shè)定的，各個系數(shù)沒有必要用同樣的量化步長。通常應(yīng)該盡量利用系統(tǒng)的硬件字長，把各個系數(shù)設(shè)置得相對量化誤差為最小。本節(jié)將先給出系數(shù)相對量化的算式，然后考慮濾波器系數(shù)的量化誤差對濾波器性能的影響。因為量化是一個非線性問題，很難有封閉形式的解，我們只能用一些數(shù)字例子說明在進(jìn)行濾波器結(jié)構(gòu)選擇時應(yīng)考慮的原則。

32系數(shù)和運算量化的影響系數(shù)相對量化的數(shù)學(xué)算式和子程序：十進(jìn)制的相對量化，就是截取有效位的長度。若xd=[3.14159,0.38455,54.3753,-0.134273]，鍵入xdr2=chop(xd,2)和xdr3=chop(xd,3)，得到的結(jié)果分別為xdr2=3.100.3854.0-0.13Xdr3=3.140.38554.4-0.134xdr2與xd之間的最大誤差為十進(jìn)制第三位的±5。所以其最大相對誤差的范圍為±(0.5~5)%。

33系數(shù)和運算量化的影響要實現(xiàn)二進(jìn)制數(shù)的相對量化，應(yīng)把原量化子程序bqtize中的步長deltax取成輸入量乘以二進(jìn)制負(fù)冪，故它成為隨各個輸入量不同的數(shù)組。得到子程序brqtize.m

functiony=brqtize(x,N)ax=abs(x);%去掉符號m=ceil(log(ax)/log(2));%ax=2^mdeltax=2.^(m-N);%求出量化步長數(shù)組xint=round(ax./deltax);%量化整數(shù)y=sign(x).*xint.*deltax;%量化輸出

34系數(shù)和運算量化的影響程序中采用了

(9.3.1)來確定數(shù)x的最高二進(jìn)制位數(shù)m，其中表示向上取整。得到m，就可以知道該數(shù)的二進(jìn)制上限，從而確定最小量化單位Δ=2(m-N)，參看式

(1.1.2)就可以理解此子程序中后兩條語句。

35系數(shù)和運算量化的影響這個子程序的優(yōu)點一是適用于x為數(shù)組的情況，不用一個一個系數(shù)去量化，特別方便。二是對所有的系數(shù)，其相對量化誤差都以各個系數(shù)本身的絕對值來計算，最大相對誤差都是2^-N。因此它可以對濾波器系數(shù)向量實現(xiàn)統(tǒng)一的量化。仍取前面給出的xd，將它量化至4位二進(jìn)制。鍵入xdq=bqtize(xd,4)得到

xdq=3.250.375056.0-0.1406求其相對誤差(xdq-xd)./abs(xd)，得到ans=0.0345-0.02480.0299-0.0473相對誤差比較均勻，都在2^(-4)=0.0625以下。

36系數(shù)和運算量化的影響例9.3.2設(shè)一個四階FIR濾波器，其傳遞函數(shù)為：算出它的級聯(lián)結(jié)構(gòu)，比較在兩種結(jié)構(gòu)下系數(shù)量化為四位二進(jìn)制對頻率響應(yīng)的影響。解：先用[sos,G]=tf2sos(b,1)求其級聯(lián)結(jié)構(gòu)參數(shù)，分別求b,sos,G的量化值bq=brqtize(b,4);sosq=brqtize(sos,4),Gq=brqtize(G,4)得出的bq=0.07030.37500.3750…然后用這幾個參數(shù)計算濾波器的幅特性，并與量化前參數(shù)產(chǎn)生的幅特性進(jìn)行比較，得出如下曲線。37系數(shù)和運算量化的影響以理想情況作為標(biāo)準(zhǔn)，可以看出，在本例中，用級聯(lián)結(jié)構(gòu)時，系數(shù)量化對系統(tǒng)幅頻特性的影響較直接結(jié)構(gòu)為小。

38系數(shù)和運算量化的影響系數(shù)量化對IIR濾波器的影響舉例IIR濾波器分母系數(shù)存在著遞歸計算的問題，它的系數(shù)誤差可能在循環(huán)計算中不斷擴(kuò)散，因此，一般地說，對系數(shù)誤差的敏感程度比FIR濾波器嚴(yán)重。特別是當(dāng)IIR濾波器的極點靠近z平面上單位圓時，幅頻特性通常要出現(xiàn)很大的峰值。系數(shù)的量化誤差造成的影響很難用解析公式表示，只能通過實例來計算和觀察。39系數(shù)和運算量化的算例9.3.2一IIR濾波器具有如下傳遞函數(shù)將系數(shù)量化為七位二進(jìn)制，分析其極點的變化。解：系數(shù)量化前、后的極點由下列核心語句求得：aq=brqtize(a,7)%量化后分母系數(shù)向量

p=roots(a),%量化前的極點pq=roots(aq),%量化后的極點其量化前、后零極點的位置畫在z平面上如下圖。鍵入abs(pq),得[1.0079,1.0079,0.9613,0.5000]40系數(shù)和運算量化的影響 量化前的零、極點用‘o’,’x’表示，而量化后的零、極點用‘+’,’*’表示。從計算的數(shù)據(jù)看出，量化后極點的模大于1，出了單位圓，系統(tǒng)由穩(wěn)定變?yōu)椴环€(wěn)定。41系數(shù)和運算量化的影響從這個例子看到，濾波器的級聯(lián)結(jié)構(gòu)對系數(shù)量化的敏感性較之直接結(jié)構(gòu)要小。高次多項式系數(shù)的微小變化往往會引起根的顯著變動。因此，階次高的濾波器不宜采用直接結(jié)構(gòu)。系數(shù)量化還會造成其它一些問題，如引起極限環(huán)振蕩等，這些問題一般出現(xiàn)在量化位數(shù)較小的單片式嵌入數(shù)字系統(tǒng)中。當(dāng)位數(shù)達(dá)到16位時，這類問題已經(jīng)不很突出了。研究系數(shù)量化問題還沒有完美的理論，主要靠針對具體問題的仿真計算來判別量化的影響。

42系數(shù)和運算量化的影響運算量化的影響運算精度的規(guī)律如下：兩個數(shù)相加，輸出的有效數(shù)位與輸入大數(shù)的有效數(shù)位同；兩數(shù)相減，輸出的相對誤差將增大，增大的倍數(shù)等于輸入大數(shù)與差之比；兩數(shù)相乘，輸出的有效數(shù)位等于輸入兩數(shù)的最低有效數(shù)位。所以只要避免兩個數(shù)量相近的大數(shù)的相減，運算器的量化誤差與輸入數(shù)據(jù)的量化誤差數(shù)量級相同。信號運算的細(xì)致過程可以用信號流圖來表示。加法器和乘法器如果工作正常（沒有溢出），應(yīng)該能夠保持運算的量化精度。43系數(shù)和運算量化的影響加法和乘法器的問題一般不在于量化，而在于發(fā)生溢出。溢出后如果不加特殊處理，加法器的輸出會從正的最大值跳變到負(fù)的最大值，或反之。這種情況有時會引起大范圍振蕩，徹底破壞濾波器的工作，那時，上面的量化模型全都無效了，所以這是絕不容許發(fā)生的。要避免這種現(xiàn)象發(fā)生，除了對輸入量的動態(tài)范圍進(jìn)行控制之外，應(yīng)該把加法器設(shè)置成溢出后飽和的工作狀態(tài)，不讓它在溢出后出現(xiàn)正負(fù)號的突跳。449.4在雙音撥號中的應(yīng)用雙音多頻(DTMF—DualToneMultiFrequency)是用按鍵進(jìn)行電話撥號的體制。它不單用在電話中，還可以用于傳輸十進(jìn)制數(shù)據(jù)的其它通信系統(tǒng)中。DTMF也廣泛應(yīng)用于電子郵件和銀行系統(tǒng)。第一章1.2節(jié)中已經(jīng)對它的功能進(jìn)行了描述，4個低頻頻率表示四行，3個高頻頻率表示三列，兩者的組合共可提供4×3＝12個字符。包括0到9的十個數(shù)字，加上字符*和#。本節(jié)要討論如何解決雙音變?yōu)閿?shù)字的問題。45在雙音撥號中的應(yīng)用

DTMF通信系統(tǒng)是一個很典型的小型信號處理系統(tǒng)，它既有模擬信號的生成和傳輸部分（這要用到D/A轉(zhuǎn)換）；又有把它轉(zhuǎn)為數(shù)字信號（這要用到A/D轉(zhuǎn)換）并進(jìn)行數(shù)字處理的部分；而為了提高系統(tǒng)的檢測速度和降低成本，還開發(fā)了一種特殊的DFT算法，稱為Goertzel算法；這種算法在國外的幾乎每一本數(shù)字信號處理的教材上都要介紹，說明了人們對它的重視。這種算法既可以用硬件（專用芯片）、也可以用軟件實現(xiàn)，所以DTMF系統(tǒng)的設(shè)計問題是理論與工程相結(jié)合的一個很好的典范。46在雙音撥號中的應(yīng)用先研究雙音信號的生成問題：DTMF的兩個音頻可以用計算法或查表法產(chǎn)生。用計算法得到正弦波形的缺點是要占一些運算時間；查表法的速度較快，缺點是要占一定的存貯空間；兩個正弦波的數(shù)字樣本按比例相加在一起。因為采樣頻率是8KHz，硬件必須每125毫秒輸出一個樣本。將這個疊合信號送到D/A變換器變換成模擬音頻信號，通過電話線路傳送到交換機(jī)。

47在雙音撥號中的應(yīng)用在接收端，將收到的模擬音頻信號進(jìn)行A/D變換，恢復(fù)為數(shù)字信號，然后檢測其中的音頻頻譜來確定所發(fā)送的數(shù)字。檢測算法可以用FFT，也可用一組濾波器來提取所需頻率。當(dāng)要檢測的音頻數(shù)目比較少時，用濾波器組更節(jié)省硬件。如果用FFT算法實現(xiàn)該DFT計算，計算量(復(fù)數(shù)乘法和加法)是Nlog2N。好處是立即得到DFT的所有N個值，但至少要N個存儲器。然而，如果只希望計算DFT的K個點,而K<<N，則直接計算可以節(jié)省很多內(nèi)存。下面介紹的Goertzel算法就屬于后者。48DFT的Goertzel算法Goertzel算法利用相位因子Wnk的周期性，將DFT運算表示為線性濾波運算，由于WN-kN=1，可用該因子去乘DFT，則它是兩個序列的卷積。一是長度為N的輸入序列x(n)，另一則是脈沖響應(yīng)序列 可以看作x(n)通過具有該脈沖響應(yīng)的濾波器。該濾波器在n＝N點的輸出就是頻點k處的DFT.49DFT的Goertzel算法單位脈沖響應(yīng)為hk(n)的濾波器的系統(tǒng)函數(shù)就是(9.4.4)式的z變換，容易求出為這個濾波器只有一個位于單位圓上頻率為ωk=2πk/N處的極點。Goertzel算法的好處不在于節(jié)省時間，而在于節(jié)省空間。如果只需要K個DFT樣本，可以只用K個并行的單極點濾波器來分別計算這K個樣本，，這就可以大大節(jié)省硬軟件資源。50DFT的Goertzel算法根據(jù)這個濾波器的差分方程，因為可以用迭代方法計算yk(n)，如下圖。要執(zhí)行該計算，可以只算一次相位因子WN-k，將其存貯起來。迭代N次后，預(yù)期的輸出為：X(k)=yk(N)。51DFT的Goertzel算法上式中包含對 的復(fù)數(shù)運算，想避免它，可將成對的復(fù)共軛極點 組合在一起計算，導(dǎo)出雙極點濾波器的系統(tǒng)函數(shù):它可以看成用兩個差分方程組的聯(lián)立:

初始條件為

vk(-1)=vk(-2)=0。運算結(jié)構(gòu)圖見后：52DFT的Goertzel算法53DFT的Goertzel算法(9.4.9)中的遞推關(guān)系對n=0，1，……，N重復(fù)N+1次，每次只需要計算一次實數(shù)乘和兩次實數(shù)加，而帶有復(fù)數(shù)運算的方程(9.4.10)僅在n＝N時刻計算一次。所以，對實數(shù)序列x(n)，由于對稱性，用這種算法求出X(k)和X(N-k)的值只需要N次實數(shù)乘法和一次復(fù)數(shù)乘法運算。實際上，不必計算復(fù)數(shù)值X(k)，只要求出幅度平方值|X(k)|2，故計算的最后一步還可以簡化，從而完全避開復(fù)數(shù)運算。54Goertzel算法子程序按(9.4.7)式編成子程序gfft。它根據(jù)輸入序列x，和規(guī)定的DFT樣本序號k，計算DFT樣本X。functionX=gfft(x,k)%用Goertzel算法計算序號為k的DFT樣本N=length(x);x1=[x,0];%遞推要N+1次，d1=2*cos(2*pi*k/N);%濾波中間項系數(shù)v=filter(1,[1,-d1,1],x1);%用濾波實現(xiàn)卷積W=exp(-i*2*pi*k/N);%為下一步求WX=v(N+1)-W*v(N);%求出第k個DFT55在雙音撥號中的應(yīng)用DTMF信號參數(shù)選擇的考慮因素：(1).為了抗干擾，要檢測8個頻率及其倍頻；(2).為了能分辨清這16個頻率，要求采集數(shù)據(jù)的長度N>Fs/D，已知Fs=8000，而根據(jù)最小相鄰基頻間隔（73Hz）大于2D的要求，D<37，得到N應(yīng)在200以上。而為了提高檢測速度，又希望N盡量小些。(3).為了檢測盡量精確，希望這16個離散頻點靠近k=整數(shù)處，這需要仔細(xì)選擇N。最后選定的規(guī)范值是N=205。56在雙音撥號中的應(yīng)用本節(jié)用一個MATLAB程序hc941來演示DTMF雙頻撥號的全過程。包括以下幾部分：（1）首先程序要求用戶輸入一個數(shù)字，然后根據(jù)這個數(shù)字查出它對應(yīng)的兩個頻率，并生成相應(yīng)的雙頻信號x(n)，并產(chǎn)生聲音作為標(biāo)志。（2）在接收端，程序截取收到的雙頻信號并進(jìn)行傅立葉變換，求出它在八個規(guī)定頻率上的樣本幅度，畫出頻譜圖。（3）從頻譜中分別取出滿足規(guī)定電平的高低頻分量的兩個下標(biāo)。再由下標(biāo)找到它的ASCII碼和數(shù)字。57程序hc941的說明此程序分為三段。第一段是雙頻信號的生成，從輸入一個數(shù)字開始，以產(chǎn)生相應(yīng)的雙頻信號結(jié)尾；第二段是接收端信號處理，從收到雙頻信號開始，進(jìn)行截斷和量化，并計算八個基頻DFT，以畫出DFT幅度分布圖結(jié)尾；第三段是根據(jù)基頻幅度分布圖,截取其兩個峰值，查找和顯示輸入字符。各段之間用pause命令隔開。58d=input(‘鍵入一位電話號碼=’,‘s’);symbol=abs(d);%求它的ASCII碼tm=[49,50,51,65;52,53,54,66;55,56,57,67;42,48,35,68];%16個鍵的ASCII碼，檢測輸入碼與誰相符forp=1:4;forq=1:4;iftm(p,q)==abs(d);break,end%檢測列號qendiftm(p,q)==abs(d);break,end%檢測行號pend程序hc941的說明59f1=[697,770,852,941];%行頻率向量f2=[1209,1336,1477,1633];%列頻率向量n=0:2040;%為了發(fā)聲，加長序列%構(gòu)成雙頻信號x=sin(2*pi*n*f1(p)/8000)+sin(2*pi*n*f2(q)/8000);sound(x);%發(fā)出聲音disp('雙頻信號已經(jīng)生成并發(fā)出'),pause％以上是第一段程序，表示發(fā)送端的功能。程序hc941的說明60程序hc941的說明（第二段）接收檢測端的計算8點DFT的程序N=205; %截取信號的長度 k=[1820222431343842];%待求的DFT序號form=1:8;X(m)=gfft(x(1:N),k(m));%算八點DFT值endval=abs(X);stem(k,val,'.');%畫出八點DFT向量disp('圖上顯示檢測到的八個基頻的DFT幅度');

％信號處理工具箱提供的goertzel函數(shù)能夠取代帶下劃線的三行程序。61程序hc941的說明62在雙音撥號中的應(yīng)用（第三段程序，選出兩個峰值基頻，查數(shù)字）limit=80;%設(shè)定門限，這個數(shù)與N有關(guān)fors=5:8;ifval(s)>limit,break,end%查找列號endforr=1:4;ifval(r)>limit,break,end%查找行號enddisp(['接收端檢測到的號碼為',setstr(tm(r,s-4))])%顯示接收到的字符63在雙音撥號中的應(yīng)用在信號處理工具箱版本5.1中，也給出了Goertzel算法的子程序goertzel.m。它的調(diào)用格式與gfft相似，不同在于第二變元k可以是數(shù)組。所以編程時可以省去一個for循環(huán)。另外它的k指的是fft的下標(biāo)，所以它等于頻率的fft編號加一。因此在程序hc941中，有波紋線的三行for循環(huán)語句可以換成語句X=goertzel(x(1:N),k+1)。第三段程序中門限值limit的設(shè)定與Goertzel算法所得的|H(k)|應(yīng)得值有關(guān)，理想值是N/2，實際應(yīng)略小些。若取N不同，此門限也必須修改。649.5正余弦信號的譜分析

數(shù)字信號處理的一個重要用途是在離散時間域中確定一個連續(xù)時間信號的頻譜，通常稱為頻譜分析。如果連續(xù)時間信號ga(t)是限帶的，那末它的離散時間等效物g(n)應(yīng)當(dāng)能給出ga(t)頻譜的近似的估計。然而，這種近似取決于g(n)在多大程度上與ga(t)相近，在大多數(shù)情況下，ga(t)是在-∞<t<∞范圍內(nèi)定義的，要估計一個無限長信號的頻譜是不可能的，g(n)必須有有限的數(shù)據(jù)長度，且定義在有限范圍內(nèi)。這就要涉及很多因素：如采樣的周期T，截取信號的長度N等。65正余弦信號的譜分析假定表征正余弦信號的基本參數(shù)，如振幅、頻率和相位，不隨時間故變，則此信號g(n)的付立葉變換G(ejω)可以用計算它的DTFT得到

實際上，把g(n)首先乘以一個長度為N的窗函數(shù)w(n)，使它變成一個長為N的有限序列，g1(n)=g(n)w(n)，對g1(n)求出的DTFTG1(ejω)，把它作為原連續(xù)模擬信號ga(t)的頻譜近似估計。66正余弦信號的譜分析為了快速計算，要用DFT來求DTFT。故在0≤ω≤2π區(qū)間等分為R點，求出離散付立葉變換DFT。為了保證足夠的分辨率，DFT的長度R選得比窗長度N大，其方法是在截斷了的序列后面補(bǔ)上R-N個零。計算采用FFT算法。我們更詳細(xì)地考察一下上面的方法。特別要分析加窗的效果，以及和由DFT樣本來估計DTFT頻率采樣值的問題。

67正余弦信號的譜分析在討論由G1(k)來估計頻譜G1(ejω)和G(ejω)時，需要探討一下這些變換和它們所對應(yīng)的頻率之間的關(guān)系，R點的DFTG1(k)與它的DTFTG1(ejω)的關(guān)系為：

0≤k≤R-1 (9.5.2)歸一化的數(shù)字頻率ωk和DFT樣本序號k的關(guān)系為

(9.5.3)模擬角頻率和k的關(guān)系為(其中T是采樣周期) (9.5.4)68正余弦信號的譜分析首先考慮單頻率正余弦序列的頻域分析。設(shè)一個具有數(shù)字角頻率ω0的余弦信號為：

(9.5.5)把這個序列表為指數(shù)序列：

(9.5.6)查表得知它的DTFT為：

因此，G(ejω)是一個以2π為周期的ω的周期信號，每個周期中在±ω0處有兩個沖擊信號。69正余弦信號的譜分析取g(n)的一個有限長序列假如g1是一個頻率為10Hz而采樣頻率為64Hz的32點序列，則用計算其DFT的程序如下：

N=32;f=10;Fs=64l;n=0:N-1; g1=cos(2*pi*f*n/Fs);k=n; G1=fft(g1,N);stem(k,abs(G1),’.’)%畫出DFT G=fft(g1,1024);plot(N*…)%畫出DTFT曲線計算結(jié)果見圖9.5.1上圖。它的DFT只有兩個點不等于零，位于k=5和k=27處。

70正余弦信號的譜分析71正余弦信號的譜分析在上圖中，k=5對應(yīng)于頻率10Hz，

k=27對應(yīng)于頻率54Hz（也就是-10Hz），這樣DFT確實正確地分辨了余弦信號的頻率。但是這樣理想的結(jié)果是碰巧得到的，因為我們恰好截取了五個完整的余弦周期(f*N/Fs=5)。如果截取的不是整數(shù)周期，情況就不同了。例如把頻率f改為11Hz，而采樣頻率仍為64Hz和窗長度仍為32點，用同樣的程序計算此余弦信號的頻譜，則計算結(jié)果見圖9.5.1下圖。頻譜圖上k=5和k=27處都有較大的峰值，而其它的點上幅度不再為零。這樣的現(xiàn)象稱為頻率泄漏，來源于截斷效應(yīng)。72正余弦信號的譜分析

從圖9.5.1中可以看到，上下兩個子圖的DTFT形狀是很相像的。兩個DFT樣本點所以有那末大的差別，原因就在于上圖的采樣點位置正巧都在頻譜的零點。實際工程中，輸入信號的參數(shù)是未知的，而且通常含有豐富的各種頻譜，這種理想情況不會出現(xiàn)。為了理解圖中的頻譜形狀，可以把有限序列g(shù)1(n)看作無限序列g(shù)(n)和長度為N的矩形窗序列的乘積。乘積的頻譜應(yīng)該等于乘子頻譜的卷積。結(jié)果，序列g(shù)1(n)的DTFTG1(ejω)是窗函數(shù)wd(n)的DTFTWdR(ejω)移頻±ωO并加權(quán)后的頻譜的的和。

73正余弦信號的譜分析長度為N的矩形窗序列wd(n)的DTFT可從圖3.2.2推論得到。其DTFTWdR(ejω)

形狀如下圖所以移頻±ωO疊加就得到前面截斷正余弦信號的頻譜。74正余弦信號的譜分析對于頻率為11Hz的序列，它的數(shù)字頻率為ω0=11/64*2π=0.344π，因序列長度N=32，頻率分辨率為Δω=2π/N≈0.0625π。ω0

/Δω=5.5,因此ω0在k橫坐標(biāo)上是在5和6之間。DFT是由長度32的矩形窗的頻譜，向左右各移動5.5后再相加并乘1/2而得的。在數(shù)字角頻率0~2π(數(shù)字下標(biāo)k=0~31)范圍內(nèi)，存在著兩個尖峰，一個在k=5.5處，表現(xiàn)為k=5和k=6之間兩處的DFT，另一個在k=32-5.5= 26.5處，表現(xiàn)為k=26和k=27處的兩個DFT。所有其它的DFT樣本是由窗函數(shù)的泄漏引起的DTFT的旁瓣所造成。75正余弦信號譜分析算例9.5.2當(dāng)輸入信號有一個以上的正余弦分量時，上述問題變得更加嚴(yán)重。例：9.5.2考察DFT的長度對雙頻率信號頻譜分析的影響。設(shè)待分析的信號為X(n)=0.5sin(2πf1n)+sin(2πf2n)0≦n≦N-1令兩個長度為N=16的正余弦序列的數(shù)字頻率為f1=0.22及f2=0.34。解：現(xiàn)在計算x(n)在取不同R值時的DFT值。為此采用下列程序，其中R,N,f1,f2都是可設(shè)定的。76正余弦信號譜分析算例9.5.2程序hc952N=input('信號長度');%R=input('DFT長度');%fr=input('[f1,f2]=');%輸入兩個頻率n=0:N-1;%設(shè)定自變量向量x=0.5*sin(2*pi*n*fr(1))+sin(2*pi*n*fr(2));

%兩個正弦信號合成X=fft(x,R);k=0:R-1;%求其DFTstem(k,abs(X),'.');gridon%畫圖77正余弦信號譜分析算例9.5.2取R為四個不同值16,32,64,128，畫出的四個DFT樣本幅度圖|X(k)|如圖9.5.2上面四個子圖。從第一個子圖上，很難看出它有兩個峰值。因此就要提高它的分辨率。把DFT的長度R由16增加到32，64和128,逐漸可以看出它有兩個峰值。注意這幾個圖的橫坐標(biāo)為k，可換算為數(shù)字頻率f=ω/2π=k/R。這樣就能確定峰值的位置大體在f=0.21和0.35附近。另外，增加DFT的長度R減小了相鄰樣本間的頻率間距，提高頻譜的視在分辨率，因而可以提高樣本位置的測定精度。78正余弦信號的譜分析79正余弦信號譜分析算例9.5.2上例中兩對沖擊函數(shù)，本來應(yīng)該可以分辨的，但由于和有限長的窗函數(shù)的頻譜進(jìn)行卷積，形成了較寬的頻譜圖形。所有的窗函數(shù)主瓣寬度都與N成反比。因此增加序列長度N可以有效地提高頻譜的實際分辨率。雖然矩形窗的主瓣寬度最小，但它的旁瓣幅度太大，造成嚴(yán)重的頻率泄漏，這會使頻譜分析的可靠性和精確性下降。因此，人們寧可選一個沒有旁瓣的窗函數(shù)，而靠加大N來提高頻譜的分辨率。下面的例子說明用Hamming窗對改善分辨率的影響。80正余弦信號譜分析算例9.5.2例9.5.3在上例中若把兩個正弦波的頻率取得較近，令fr=[0.22,0.25]，試問應(yīng)該怎樣選參數(shù)才能在頻譜分析中分辨出這兩個分量？解：要使這兩個頻率之間隔開一個樣本，分辨率至少應(yīng)達(dá)到Δf=0.03/2=0.015。因為此處的數(shù)字頻率是對采樣頻率Fs進(jìn)行歸一化的，即fr的最大值為1。因此總的樣本數(shù)N至少要達(dá)到1/0.015=66。為了分析加窗的影響，只要在程序hc952a中增加一句x1=x.*hamming(N)’，就構(gòu)成hc953的核心語句，然后對x1求DFT即可。81正余弦信號譜分析算例9.5.2程序hc953選擇不同N,R和窗函數(shù)運行后，可得到圖9.5.3。從圖上可以看出，信號的兩個峰值是能夠分辨的。將橫坐標(biāo)放大，可以得到圖9.5.4。左邊的一列是加矩形窗的結(jié)果，右邊的一列是加哈明窗的結(jié)果，可以看出，使用無旁瓣的窗函數(shù)得到的頻譜函數(shù)比較光滑，便于分辨峰值位置和準(zhǔn)確的數(shù)值。所以在頻譜分析儀中都要采取加窗的措施。此外，為了提高實際的分辨率，應(yīng)該盡量增加信號長度N及DFT長度R(≥N)，82正余弦信號的譜分析83正余弦信號的譜分析84非平穩(wěn)信號的譜分析在工程實踐中，往往不可能任意增加數(shù)據(jù)長度N。人們所遇到的大多數(shù)信號，例如語音通信、雷達(dá)信號…等等，都不是平穩(wěn)的。在這種情況下，增加取樣的長度N未必有實際的意義，很多的時候，要求人們進(jìn)行非平穩(wěn)信號的頻譜分析，于是出現(xiàn)了多種多樣的分析方法。這些都要用比較高深的數(shù)學(xué)工具，這里簡單介紹一種，稱為短時間傅立葉變換(STFT－Short-TimeFourierTransfom)方法。。85非平穩(wěn)信號的譜分析序列x(n)的STFT定義為

(9.5.12)

其中，wd(m)是一個適當(dāng)選擇的窗函數(shù)序列。這個式子說明，窗函數(shù)wd(m)在x(n)中選擇一個有限長度序列。在這段時間內(nèi)，近似把非平穩(wěn)過程當(dāng)作平穩(wěn)過程來進(jìn)行傅立葉分析。對于一個非平穩(wěn)序列，如果窗函數(shù)隨時間移動，則得出的頻譜也將隨時間而變化。對于隨時間而變化的頻譜.時間軸和頻率軸就占了兩維，表示其幅度的第三維只能用顏色或圖形的灰度來表示。86非平穩(wěn)信號譜分析算例9.5.4程序hc954可以給出一個線性調(diào)頻(Chirp)信號的STFT頻譜圖，STFT函數(shù)specgram的缺省輸入變元為，F(xiàn)FT長度R=256，序列長度N=R/2，窗函數(shù)為漢寧窗。t=0:0.001:2;%1kHz采樣頻率，持續(xù)2秒x=chirp(t,0,1,150);%1秒內(nèi)從直流達(dá)到150Hz的線性調(diào)頻信號

subplot(2,1,1),plot(t,x)%畫出信號曲線

subplot(2,1,2),specgram(x);%顯示頻譜圖運行此程序得到圖9.5.5。87非平穩(wěn)信號譜分析算例9.5.4上圖為調(diào)頻信號時間曲線。下圖中橫坐標(biāo)為時間軸，表示窗的位置；縱坐標(biāo)為頻率，向上為頻率增加；灰度表示幅特性，全白意味著幅度最大。向上的白色斜線說明峰值頻率隨時間的增加向高頻端移動。

889.6音樂信號處理由于CD、DVD和家庭視聽技術(shù)的發(fā)展，音樂的錄制和加工已經(jīng)愈來愈多地采用數(shù)字技術(shù)。音樂的錄制大致分以下幾個步驟：首先在一個隔音的舞臺上把樂隊中各個樂器的聲音分別錄在一個多磁道磁帶的各個獨立磁道上；然后由音響工程師把各個磁道上的信號進(jìn)行單獨的處理，加入特定的聲音效果；最后在一個混音系統(tǒng)中把這些信號進(jìn)行最后的合成，錄制在一個立體聲的雙磁道磁帶上。這些特定的聲音效果已愈來愈多地靠數(shù)字信號處理技術(shù)來實現(xiàn)，本節(jié)將作一簡單的介紹。89音樂信號處理時域處理方法在一個像音樂廳那樣的封閉空間中，人們聽到的聲音包括直接傳播聲音、一次反射聲音和混響(reverberation)等幾種成分。一次反射聲音主要是由直接傳播聲音在近距離內(nèi)的直接回聲構(gòu)成，而混響則是由密集的回聲形成的。所以在一個隔音的舞臺上錄制的音樂與音樂廳中錄制的效果就不一樣，人們聽起來會感到‘不自然’。這時就需要用數(shù)字濾波器來人為地改變錄制的信號，通過增加一些回聲，使它接近于音樂廳中的效果。90音樂信號處理回聲可以用延遲單元來生成。直接聲音和它的延遲了R個周期的單個回聲可以用如下的差分方程表示：

(9.6.1)其中α<1表示回聲的衰減系數(shù)。上述差分方程也可以用傳遞函數(shù)為: (9.6.2)的FIR濾波器來實現(xiàn)，這實際上是一個梳狀濾波器，其結(jié)構(gòu)如圖9.6.1(a)。91音樂信號處理可以用實驗來證明這一點。將‘大家好’文件中的變量x作為輸入語聲x(n)，讓它通過(9.6.2)的濾波器，把得到的y(n)放音檢驗，為了得到0.3秒的延遲和30%的衰減，濾波器參數(shù)R選為0.3Fs=0.3×22050=6615，α=0.3。MATLAB程序為hc961如下：loaddajiahaox1=[x;zeros(5000,1)];%把x加長y=filter([1,zeros(1,6615),0.3],1,x1);

％濾波輸出sound(y,22050) ％發(fā)聲為了使輸出時間足夠長，把x加了5000個樣本。92音樂信號處理多重回聲可以用一個IIR濾波器來實現(xiàn)。其傳遞函數(shù)為

(9.6.4)在上述程序第二行中，把濾波器的分子分母系數(shù)掉換一下，同時把α取大一些，又把x加得更長，寫成以下求合成聲音y1的語句x1=[x;zeros(50000,1)]y1=filter(1,[1,zeros(1,6615),0.8],x1);sound(y1)這樣就可以檢驗多重回聲效果。93音樂信號處理圖示為這兩種梳狀濾波器的結(jié)構(gòu)。不過它們還無法滿足音樂處理的要求。首先是因為這種濾波器的幅頻特性不是常數(shù)，對不同頻率的聲音諧波響應(yīng)不均勻；其次是這種回聲太單調(diào)，每秒中的回聲數(shù)目太少會導(dǎo)致聲音的顫動。

94音樂信號處理頻域處理方法把分別錄制的各種樂器或歌手的聲音進(jìn)行混合時，通常要由音樂工程師修改它們的頻率響應(yīng)。方法是讓信號通過一個均衡器（equalizer），其作用是使這些聲音在混合信號中充分表現(xiàn)。另外也需要通過‘?dāng)U大’或‘削減’在某些頻率范圍的信號，以修正低頻和高頻的信號之間的關(guān)系。通常用許多一階和二階的參數(shù)可調(diào)的濾波器級聯(lián)起來實現(xiàn)這個功能。濾波器的結(jié)構(gòu)選擇的主要要求之一是調(diào)整方便，最好是調(diào)一個參數(shù)只影響一個應(yīng)用指標(biāo)，而且可調(diào)參數(shù)要少。95音樂信號處理一階均衡器的可調(diào)濾波器結(jié)構(gòu)圖如圖9.6.1(a)，它到兩個輸出端的傳遞函數(shù)分別為

(9.6.6)其中96音樂信號處理若把低頻輸出乘以K與高頻輸出相加，如右圖，就成為低頻均衡濾波器；其傳遞函數(shù)為

(9.6.8)在α和K給定時,其頻率響應(yīng)可用程序hc962計算:K=2.5;alpha=0.7;b1=[alpha,-1];a1=[1,-alpha];b=polyadd(0.5*(1+K)*a1,b1);[H,w]=freqz(b,a1);設(shè)定若干個alpha值，得出的曲線如下左圖。97音樂信號處理設(shè)定的K值小于1時，得出的曲線如下右圖。因此，可以根據(jù)需要，設(shè)置可調(diào)的電位器，來提高低頻的增益或高頻的增益，并可調(diào)節(jié)邊界頻率。98音樂信號處理二階均衡器也采用圖9.6.1的結(jié)構(gòu)圖，不同點僅僅在于把其中的A1(z)換成A2(z)。

把A2(z)替換了A1(z)后，公式(9.6.6)到(9.6.9)全部適用。因此可以寫出

(9.6.11)它的頻率響應(yīng)見圖9.6.5。

99音樂信號處理二階均衡器的可調(diào)參數(shù)有三個：K,α和β。可以分別調(diào)節(jié)幅度、帶寬和中心頻率。通常不用更高階的均衡器，寧可用多個均衡器級聯(lián)。1009.7變采樣率數(shù)字濾波

前面所討論的信號處理的各種方法都是把采樣率Fs視為固定值，即在一個數(shù)字中只有一個采樣頻率。但在實際系統(tǒng)中，經(jīng)常會遇到采樣率的轉(zhuǎn)換問題，即要求一個數(shù)字系統(tǒng)能工作在“多采樣率”狀態(tài)。在數(shù)字系統(tǒng)愈來愈普及的情況下，各個數(shù)字系統(tǒng)都有自己不同的標(biāo)準(zhǔn)，它們之間的銜接也不可避免地要遇到這個問題。近年來，建立在采樣率轉(zhuǎn)換基礎(chǔ)上的“多采樣率數(shù)字信號處理”已成為數(shù)字信號處理學(xué)科中的主要內(nèi)容之一。

101變采樣率數(shù)字濾波直觀地想，變換采樣率的最簡單方法是：先將以采樣率Fs1采集的數(shù)字信號進(jìn)行D／A轉(zhuǎn)換，變成模擬信號；再按采樣率Fs2進(jìn)行A／D變換，這就實現(xiàn)了從Fs1到Fs2的采樣率轉(zhuǎn)換。但這樣做系統(tǒng)復(fù)雜，且易使信號受到額外的干擾和崎變，所以在實用的變采樣率系統(tǒng)中，改變采樣率并不經(jīng)過模擬信號，而完全是在數(shù)字域?qū)崿F(xiàn)的。

102變采樣率數(shù)字濾波采樣率轉(zhuǎn)換通常分為“抽取(Decimation)”和“插值(Interpolation)”。抽取是降低采樣率(所以在英文中也稱為Down-sampling)，以去掉多余數(shù)據(jù)樣本的過程，而插值則是提高采樣率(Up-sampling)以增加數(shù)據(jù)樣本的過程。變采樣率是一個非線性問題，它涉及相當(dāng)復(fù)雜的數(shù)學(xué)推導(dǎo)，通常不在本科階段討論它的設(shè)計問題。本書主要介紹抽取和插值的一般概念，一方面可以加深對采樣定理的理解，同時可以幫助讀者了解更多的應(yīng)用問題，開擴(kuò)思路。

103內(nèi)插過程的數(shù)字濾波信號的整數(shù)倍內(nèi)插整數(shù)倍內(nèi)插是在已知的相鄰兩個采樣點之間插入L-1個樣本點。由于這些樣本點并非已知的值，所以關(guān)鍵問題是如何求出這L-1個樣本值。前提是要直接在數(shù)字域完成這個插值。用數(shù)字方法實現(xiàn)整數(shù)內(nèi)插的步驟如下：在已知采樣序列x(n1T1)的相鄰兩個樣點之間等間隔插入L-1個零值樣本點，得到輸出y，則有104內(nèi)插過程的數(shù)字濾波示例式中L為大于1的整數(shù)，稱為插入因子。按上式y(tǒng)(n)是取值跳動的序列，讓它經(jīng)H(z)進(jìn)行低通數(shù)字濾波，得到教平滑的內(nèi)插結(jié)果。這種內(nèi)插方案的框圖如圖9.7.1所示。各點的波形和頻譜見圖9.7.2。105內(nèi)插中的數(shù)字濾波示例9.7.1Nx=10;L=3;fsx=2;T=1/fsx;alpha=0.5;nx=0:Nx-1;%設(shè)定x的自變量向量x=exp(-alpha*nx*T);%給出x序列的值Ny=Nx*L;ny=0:Ny-1;%設(shè)定y的自變量向量p=ones(1,L*Nx);%生成p序列y=zeros(1,L*Nx);%先y序列初始化（全置零）y(1:L:L*Nx)=x;%內(nèi)插給出y在x對應(yīng)點處的值b=[0,0.0019,0.0203,0.0772,0.1589,0.2,0.1589,0.0772,0.0203,0.0019,0];%用例733的濾波器系數(shù)y1=filter(b,1,y);%進(jìn)行濾波106內(nèi)插中的數(shù)字濾波示例9.7.1y2=exp(-alpha*ny*T/L)/L;%理想插值X=fftshift(fft(x));%x頻譜，移到對稱位置Y=fftshift(fft(y));%y頻譜，移到對稱位置P=fftshift(fft(p));%p頻譜，移到對稱位置Y1=fftshift(fft(y1));%y1頻譜，移到對稱位置Y2=fftshift(fft(y2));%y2頻譜，移到對稱位置nxm=floor(nx-Nx/2+0.5);%生成X頻率位置序列nym=floor(ny-Ny/2+0.5);%生成Y頻率位置序列%以下繪圖語句，得出各組曲線如下圖。107變采樣率數(shù)字濾波108內(nèi)插中的數(shù)字濾波示例9.7.1現(xiàn)在分析整個整數(shù)倍內(nèi)插過程中的信號變換關(guān)系。在由x變?yōu)閥時，波形沒有變，但采樣率提高了，這相當(dāng)于把序列x與一個單位樣本序列p相乘，該p的采樣頻率為Fsy=L×Fsx，生成的y序列采樣頻率與p相同，每隔L-1個點有一個非零脈沖，其值等于該時刻的x(n)。時域的相乘對應(yīng)于頻域的卷積，因此y的頻譜Y相當(dāng)于把原來x的頻譜X與頻譜P卷積，頻譜P是把奈奎斯特頻段平分的L個脈沖，時域的相乘對應(yīng)于頻域的卷積，因此y的頻譜Y相當(dāng)于把原來x的頻譜X與頻譜P卷積，頻譜P是把x的奈奎斯特頻段擴(kuò)展L倍的沖激頻譜。

109抽取過程的數(shù)字濾波信號的整數(shù)倍抽取整數(shù)倍抽取是在已知序列x中每隔M-1個采樣點取出一個樣本點，組成新的序列y。因此y序列中的所有值都是已知的，可以直接寫出下列表達(dá)式

也就是說，輸入序列x中下標(biāo)為M的整倍數(shù)的點都被取入y序列，而其它的

樣本都被舍去不用了。用MATLAB語言來表達(dá)，為：Ny=ceil(Nx/M);nx=0:Nx-1; y=x(0:M:Nx);ny=0:Ny-1;

110抽取過程的數(shù)字濾波 抽取的數(shù)學(xué)關(guān)系看起來較插入簡單。但由于采樣率的降低，造成的頻率泄漏可能變得相當(dāng)嚴(yán)重。要解決好這個問題，在進(jìn)行抽取之前，必須對原來的序列進(jìn)行預(yù)先數(shù)字濾波，稱為‘抗泄漏濾波’。它的框圖可表為圖9.7.3

111抽取中的數(shù)字濾波算例9.7.2與框圖對應(yīng)的程序為hc972，先用濾波器濾波得到x1，再抽取得到y(tǒng)1。故其核心語句為：nx=0:Nx-1;x=exp(-alpha*nx*T);%給出x序列b=[0,0.0019,0.0203,0.0772,0.1589,0.2,0.1589,0.0772,0.0203,0.0019,0];%例733的濾波器系數(shù)x1=filter(b,1,x);%進(jìn)行濾波，得出x1Ny=ceil(Nx/M);ny=0:Ny-1;y1=x1(1:M:Nx);%實行抽取，得出y1序列以下分別畫出x,x1,y1,y2的頻譜見下圖。112變采樣率數(shù)字濾波113抽取中的數(shù)字濾波算例9.7.2從圖中左邊一列看出抽取過程的時域波形，指數(shù)下降的輸入信號x經(jīng)過低通濾波成為波形x1，然后經(jīng)過抽取而降低了采樣率,成為y1.最下面的圖表示不經(jīng)過濾波直接抽取的結(jié)果。再從右邊的頻譜圖上可以看出濾波的必要性。直接抽取波形的頻譜abs(Y2)在奈奎斯特頻率處（也就是頻譜圖的左右邊界±π處）有很高的幅度，這意味著存在著很大的頻率泄漏。而經(jīng)過濾波再抽取的信號頻譜abs(Y1)在該處的幅度就小得多，如果用一個好一點的濾波器，結(jié)果會更加理想。

114分?jǐn)?shù)倍變采樣率數(shù)字濾波分?jǐn)?shù)倍變采樣率和相應(yīng)的MATLAB函數(shù)如果變換前后的采樣頻率的比滿足兩個整數(shù)之比L/M，則不難想象，可以通過先作L倍內(nèi)插后作M倍抽取方法來實現(xiàn)，也可以先抽取后內(nèi)插。其濾波器的設(shè)計當(dāng)然有不同的要求和方法。MATLAB信號處理工具箱提供了一個變采樣率函數(shù)resample，用它可以把變采樣率的全過程在一步中完成。這個變采樣率函數(shù)的典型調(diào)用格式為y=resample(x,L,M)。其中L為采樣率提高的倍數(shù)，M為采樣率降低的倍數(shù)，L和M必須取整數(shù)。

115變采樣率數(shù)字濾波的應(yīng)用整數(shù)倍抽取和內(nèi)插的應(yīng)用舉例在數(shù)字電話系統(tǒng)中，采樣頻率取為F＝8kHz，希望傳輸盡量接近4kHz的音頻帶寬。但送話器發(fā)出的信號x（t）的帶寬比4kHz大很多。因此，在A／D變換之前要對其進(jìn)行模擬預(yù)濾波，以防止采樣后發(fā)生頻率混疊失真。這就要求該濾波器在4kHz處幅頻特性很低，而為了使能傳送的信號頻帶盡量寬，又要求該濾波器的通帶接近于4kHz(例如3.8kHz)，給過渡帶寬留的裕度只有幾百Hz，這樣的要求對用模擬低通濾波器是很難做到的。

116變采樣率數(shù)字濾波的應(yīng)用為了降低對模擬預(yù)濾波器的技術(shù)要求，新方案先用較高的采樣率進(jìn)行采樣，比如取采樣率Fsx＝16kHz，這樣，模擬預(yù)濾波器的過渡帶可以從3.8kHz到8kHz，這就容易設(shè)計了。經(jīng)過A/D后，再經(jīng)M=2倍抽取，把采樣率降低至8kHz。此時問題的難點變成了抽取時的濾波器設(shè)計技術(shù)。因為這是是在數(shù)字域的濾波，用FIR結(jié)構(gòu)的數(shù)字濾波器不難設(shè)計成線形相位和過渡帶很窄的特性。這種方案最終可有效地利用通信帶寬，并增加信號數(shù)據(jù)傳送量。1179.8稀疏天線陣列設(shè)計

本書討論的信號處理都是以時間作為自變量的。而在本例中自變量卻是一維的空間位置。從本例可以看到，利用數(shù)學(xué)模型的相似性，時域的信號處理理論和方法可以推廣到空域的信號處理中。

可以發(fā)現(xiàn)，等間隔分布的線性相位天線陣元所產(chǎn)生的遠(yuǎn)場方向圖和FIR濾波器的頻率響應(yīng)在數(shù)學(xué)模型上非常相似。這就使我們可以根據(jù)方向圖的要求來設(shè)計天線陣列。在實際中往往需要去掉其中某些陣元，稱為稀疏陣列。118稀疏天線陣列設(shè)計圖9.8.1所示為N+1個等間隔配置的方向性相同的陣元。陣元間距為d，因此它們在基線上排列的坐標(biāo)為x(n)=n×d，0≤