《DSP技術及應用》課件第6章

第6章匯編程序設計6.1算法設計方法6.2

FIR濾波器的設計6.3

IIR濾波器設計6.4快速傅立葉變換的DSP實現(xiàn)

6.1算法設計方法

一個復雜的算法，如果不考慮實時實現(xiàn)，可以考慮調用現(xiàn)有的函數(shù)實現(xiàn)。如求y=sin(x)，可以直接調用C語言的sin函數(shù)。該函數(shù)直接采用浮點運算，其內部用泰勒級數(shù)公式實現(xiàn)計算，輸入和輸出參數(shù)也都為浮點數(shù)。調用C語言的sin函數(shù)不僅簡單、快捷，而且精度高，但運算量大，要求較多的運行時間。如果時間不夠，可考慮用以下方法改進。6.1.1查表法

如果輸入變量與輸出變量可以用簡單的映射關系來描述，就可以用查表法實現(xiàn)相關函數(shù)。如sin函數(shù)(這里泛指任意一個非線性函數(shù))，若在它一個周期內先均勻采樣N個點，制作成一個表，然后將輸入變量x變換成表的偏移地址，即可實現(xiàn)，如圖6.1所示。圖6.1

sin函數(shù)一個周期的采樣值(N=16)表的偏移地址的計算方法是，將x=［0~2π)弧度變換到0~N－1之間的整數(shù)index。

index=floor(x′)

(6.2)則y=sin(x)近似等價為

y=sin_table(index)(6.3)查表法速度快，但表中兩點之間的值需用靠近它們的值代替，故精度不高。要提高精度，有多種方法。如果存儲單元數(shù)富有余地，可采用加大表的長度的方法提高精度;如果存儲單元數(shù)不足，可采用線性插值的方法提高精度。設采用線性插值計算x′和index，如圖6.2所示。圖6.2線性插值示意圖則插值點x對應的y值近似等價于

y≈y′=y(index)+［y(index+1)-y(index)］(x′-index)(6.4)式中，y(index+1)-y(index)可單獨制作一張表，而不必每次都查表計算。采用更高階插值的方法可以進一步提高數(shù)據(jù)精度，但運算量也隨之增加。6.1.2數(shù)學變換法

C54x的直接運算只有加、減、乘、移位和邏輯運算等，要實現(xiàn)其他運算，需采用數(shù)學變換將其他運算變換成

這類運算才可以實現(xiàn)，最常見的數(shù)學變換為級數(shù)展開法。如正弦用級數(shù)展開，取前5項為(6.5)式(6.5)具有迭代性。若系數(shù)用Q15格式表示，從最里面的括號開始計算，均為乘加(減)運算，可用C54x的匯編專用指令快速實現(xiàn)。余弦等其他函數(shù)也可用類似的方法處理。

正弦函數(shù)也可用迭代公式(6.6)計算，但不推薦使用，因為它容易產(chǎn)生累加誤差。

sin(nθ)=2cos(θ)sin［(n－1)θ］－sin［(n－2)θ］(6.6)例6.1用查表法產(chǎn)生線性調頻信號。

線性調頻信號(LinearFrequencyModulation，LFM)的理論公式為(6.7)其中式中,K為調頻斜率，f1為起始頻率，f2為結束頻率，T為持續(xù)時間。若用數(shù)字系統(tǒng)實現(xiàn)調頻信號，需對信號進行離散化處理。設采樣頻率為fs=1/Ts，Ts為采樣周期，則t=nTs，n=0,1,2,…,N－1，N=T/Ts為一個LSF信號周期的采樣總點數(shù)。另設一個正弦周期的離散相位數(shù)為M，離散化式(6.7)，得(6.9)式中，u(n)應是u(nTs)，因Ts為常數(shù)，書寫習慣上一般省略。公式雖然簡單，但運算量卻不小。要在一個采樣周期內算出一個采樣值u(n)，需要足夠快的CPU。一般將運算公式變換成運算量小的加減運算，盡量避免運算量大的乘除運算。

令其中(6.10)(6.11)(6.12)當n>0時，f(n)可寫成遞推公式的形式

f(n)=f(n－1)+2Kcn－Kc+Kp

(6.13)

用兩個累加器實現(xiàn)：(6.14)其中φ(0)=0，φ(0)=0。從式(6.7)可以看出，u(n)仍為正弦信號，只要按式(6.13)求出n時刻的f(n)，查正弦表即可求出u(n)(直接計算正弦值的運算量很大，一般不采用)。

如f1=450Hz，f2=150Hz，T=5s，設采樣頻率fs=8kHz，

一個正弦周期的離散相位數(shù)M=1024,可算出Kc=－4.8×10－4,Kp=57.6。

6.2

FIR濾波器的設計

數(shù)字濾波器可分為IIR濾波器和FIR濾波器。從設計理

論來講，設計數(shù)字濾波器是比較復雜的過程，但借助于

MATLAB等現(xiàn)代設計軟件，可使問題大大簡化。如果利用這些工具設計出來的結果不能滿足實踐要求，則需要根據(jù)數(shù)字濾波器的設計理論進行計算，再轉化為DSP來實現(xiàn)。

FIR濾波器的系統(tǒng)函數(shù)為(6.15)對z作反變換，可得其時域輸入輸出關系為(6.16)這是一個典型的乘加運算，可采用直接型結構實現(xiàn)，如圖6.3所示。圖6.3

FIR濾波器的直接型實現(xiàn)

FIR的實現(xiàn)結構還有級聯(lián)型、頻率抽樣型、快速卷積結構和線性相位FIR濾波器(也稱系數(shù)對稱FIR濾波器)結構等多種

形式。如果系數(shù)具有對稱性，可采用線性相位FIR濾波器結構，以節(jié)省運算時間。在DSPLIB庫函數(shù)中，一般采用級聯(lián)型結構，因為級聯(lián)型結構更容易實現(xiàn)代碼的重復使用。

要實現(xiàn)FIR濾波器，需合理使用存儲器，以獲得較高的運算效率。6.2.1直接型FIR濾波器實現(xiàn)

直接型FIR濾波器延遲單元(z－1)的實現(xiàn)不同，對存儲器的要求也不同，主要有線性緩沖區(qū)法和循環(huán)緩沖區(qū)法，用來實現(xiàn)數(shù)據(jù)延遲和更新。

1.線性緩沖區(qū)法

如果在數(shù)據(jù)存儲器定義的緩沖區(qū)是雙訪問隨機存儲器(DRAM)，在一個周期內既能讀一個操作數(shù)又能寫一個操作數(shù)，則可用線性緩沖區(qū)法實現(xiàn)。

在程序存儲器中定義濾波系數(shù)h向量，首地址設為pmad。在數(shù)據(jù)存儲器定義一個緩沖區(qū)用于存放x序列值，最新的數(shù)據(jù)存放在緩沖區(qū)頂部，最“老”的數(shù)據(jù)放在底部，將輔助寄存器指針ARx初始化并指向底部，如圖6.4所示。圖6.4線性緩沖區(qū)法的數(shù)據(jù)存儲設源累加器src初始化為0，pmad指針每計算一次就自動加1，而ARx指針初始化為緩沖區(qū)底部，每計算一次就自動減1，且數(shù)據(jù)緩沖區(qū)的數(shù)據(jù)每計算一次就自動移到下一個地址(高地址)的單元上，其運算如下。第1次運算：src=src+h(N-1)×x［n-(N-1)］，x［n-(N-1)］移出到緩沖區(qū)外；

第2次運算：src=src+h(N-2)×x［n-(N-2)］，x［n-(N-2)］移到原x［n-(N-1)］單元；

…

第N次運算：src=src+h(0)×x(n)，x(n)移到原x(n－1)單元；接下來輸入x新的采樣數(shù)據(jù)x(n+1)，然后重復上述計算。一次濾波運算需N次乘加運算。C54x有專門的硬件結構和MACD指令來實現(xiàn)算法。

MACD指令格式和執(zhí)行內容如下：

MACDSmem，pmad，src

；src=src+Smem*pmad,(Smem+1)=Smem,T=Smem其中，pmad絕對尋址在重復執(zhí)行時會自動加1并指向下一個單元。例6.2

FIR線性緩沖區(qū)法濾波子函數(shù)示例。

_fir:

PSHMST1;壓棧保護ST1

SSBXSXM

SSBXFRCT

SSBXOVM;有符號小數(shù)運算,Q15狀態(tài)

;設當前采樣值x(n)已載入到A中，更新x(n)

STM#data_buffer,AR2;data_buffer為數(shù)據(jù)緩沖區(qū)x的首地址

STLA,*AR2;更新x(n)

STM#(data_buffer+N-1),AR2;從高地址開始計算

RPTZ A,#N-1

MACD*AR2-,COFF,A;COFF為系數(shù)向量h的首地址SFTA A,-16;僅輸出高16位，函數(shù)返回值必須存放在A的低16位

RET

2.循環(huán)緩沖區(qū)法

如果在數(shù)據(jù)存儲器定義的緩沖區(qū)是單訪問隨機存儲器(SARAM)，在一個周期內就只能讀一個操作數(shù)或寫一個操作數(shù)，那么線性緩沖區(qū)法就不能實現(xiàn)，可用循環(huán)緩沖區(qū)法實現(xiàn)。如圖6.5(a)、(b)所示，在數(shù)據(jù)存儲器開辟2個存儲塊存放濾波系數(shù)向量h和輸入向量x，設h首地址(低地址，以下同)為h，用ARy指針指示系數(shù)表位置，x首地址為x，用ARx指示緩

沖區(qū)位置。在濾波運算過程中,x緩沖區(qū)的數(shù)據(jù)位置不變，僅用最新的數(shù)據(jù)更新最老的數(shù)據(jù)，如圖6.5(c)所示，用x(n+1)更新x［n－(N－1)］，降低了對存儲器的要求。圖6.5循環(huán)緩沖區(qū)示意圖運算過程如下所述。

初始化：設置循環(huán)緩沖大小寄存器BK=R(緩沖區(qū)大小)，AR0=1(循環(huán)尋址地址增量)，ARy用系數(shù)向量的首地址(h)初始化，ARx用輸入數(shù)據(jù)的首地址(x)初始化，且h和x的地址必須

滿足循環(huán)尋址要求，即它們低N位地址為0，且2N>R，運算之前源累加器初始化為0。第1次運算：src=src+h(0)×x(n)，運算后指針未超邊界，ARy→h+1,ARx→x+1；

第2次運算：src=src+h(1)×x(n-1)，運算后指針未超邊界，ARy→h+2,ARx→x+2；

…

第N次運算：src=src+h(N-1)×x［n-(N-1)］，運算后指針超邊界，ARy→h,ARx→x。

接下來輸入x新的采樣數(shù)據(jù)x(n+1)更新x［n-(N-1)］。由于數(shù)據(jù)指針ARx→x(即當前x(n)所在單元)，因此指針需后退一步才可以更新x［n-(N-1)］，然后重復上述計算。

一次濾波運算需N次乘加運算。C54x有專門的硬件結構和指令來實現(xiàn)MAC算法。MAC指令格式和執(zhí)行內容如下：MAC *AR2+0%,*AR3+0%,A例6.3

FIR循環(huán)緩沖區(qū)法濾波程序示例。

.globalstart,_FIR

.mmregs

N.SET 6;

.bss

h,N,1,8；定義一個8字的邊界,8為大于6的最小整數(shù)次冪

.bssx,N,1,8

.bssy,1

.text

start:

STM#1,AR0

STM#h,AR2

STM#x,AR3

STM#y,AR4

STM#N,BK

CALL_FIR;

B

start;

_FIR:

PSHMST1

SSBXSXM

RSBXFRCT

SSBXOVM

;當前采樣值在AL中

STL

A,*AR3;更新采樣值

RPTZA,N－1

MAC*AR2+0%,*AR3+0%,A;雙間接尋址只能使用AR2~AR5

SFTA

A,－16;若為函數(shù)輸出，只能放在AL

STL

A,*AR4

MAR*AR3-0%;指針調整

POPM

ST1

RET

.end6.2.2系數(shù)對稱FIR濾波器實現(xiàn)

線性相位FIR濾波器的系數(shù)是對稱的，即h(n)=±h(N－1－n)，n=0,1,…,N－1，利用對稱性可以進一步提高運算速度。下面以N=8的偶對稱為例來說明這類FIR濾波器在C54x

上的實現(xiàn)。展開其輸出方程，可得

y(n)=h(0)×［x(n)+x(n－7)］+h(1)×［x(n－1)+x(n－6)］

+h(2)×［x(n－2)+x(n－5)］+h(3)×［x(n－3)+x(n－4)］

C54x有一條專用指令(FIRS指令)來完成系數(shù)對稱FIR濾波器計算。格式為

FIRSXmem,Ymem,pmad;B+A(32~16)×pmad→B,((Xmem)+(Ymem))<<16→A

;如果重復計算，pmad會自動加1，指令變成單字指令

FIRS可以高效地完成系數(shù)對稱的FIR濾波運算，但由于指令流水線是先讀后寫，執(zhí)行B+A(32~16)×pmad→B時對累加器A為讀，其中A(32~16)為原來A中的值，而執(zhí)行((Xmem)+(Ymem))<<16→A時對累加器A為寫，故有一個時間錯位，本次計算的A不能在本次B+A(32~16)×pmad→B中使用。如圖6.6所示，在程序存儲器安排系數(shù)表，首地址為pmad。在數(shù)據(jù)存儲器開辟2個緩沖區(qū)，大小均為N/2=4個單元，首地址分別定義為x_new和x_old，均需滿足循環(huán)尋址要求。這里N/2=4，存儲器邊界應取8,數(shù)據(jù)的安排如圖6.6所示。圖6.6

FIRS數(shù)據(jù)存儲示意圖運算過程如下所述。

初始化：AR2用x_new初始化,AR3用x_old初始化,AR0=1(地址增量);BK=N/2(循環(huán)緩沖區(qū)的大小)；

第1次：先計算x(n)+x(n-7)并存儲到AH，AR2和AR3循環(huán)尋址，地址加1;

第2次：執(zhí)行FIRS指令，B=B+［x(n)+x(n－7)］×h(0)，x(n－1)+x(n－6)→AH；AR2和AR3同上；第3~5次：與第2次類似，最后把B的結果輸出;

第6次：調整AR2和AR3指針，將new緩沖區(qū)最老的數(shù)

據(jù)復制給old緩沖區(qū)最老的數(shù)據(jù)單元,本次將x(n－3)復制到

x(n－7)；

第7次：將最新的數(shù)據(jù)(本次為x(n+1))載入new緩沖區(qū)最老的數(shù)據(jù)單元(本次為x(n－3))并調整指針，為下次FIR濾波作準備，如圖6.7所示。圖6.7更新后的FIRS數(shù)據(jù)存儲示意圖例6.4系數(shù)對稱的FIR濾波程序示例。

.globalstart,_FIR,_input

.mmregs

.data

COFF:

.word

1,2,3,4

N.SET4;

.bssx_new,8,1,8;邊界必須比N大

.bssx_old,8,1,8

.bss_input,1;數(shù)據(jù)輸入

.bssy,1

.text

start:

STM#1,AR0

STM#x_new,AR2

STM#x_old,AR3

STM#y,AR4

STM#N,BK

_FIRS:

PSHMST1

SSBXSXM

RSBXFRCT

SSBXOVM

;當前采樣值在變量_input中

LD*(_input),A

STLA,*AR2;更新采樣值x(n)

ADD*AR2+0%,*AR3+0%,A;AH=x(n)+x(n-7)

RPTZB,#N-1

FIRS*AR2+0%,*AR3+0%,COFF

STHB,*AR4

MAR*+AR2(-2)%;指針調整

MAR*+AR3(-1)%

MVDD*AR2,*AR3+0%

;用new最新的數(shù)據(jù)更新old最老的數(shù)據(jù)

POPMST1

BFIRS;

.end6.2.3

FIR濾波系數(shù)的MATLAB輔助設計

利用MATLAB函數(shù)或圖形界面fdatool可以方便地計算FIR濾波系數(shù)。

(1)函數(shù)fir1：用窗函數(shù)設計法計算FIR濾波系數(shù)，調用格式有以下三種：

B=fir1(N,wc)

B=fir1(N,wc,′highorstop′)

B=fir1(N,wc,win)

B=fir1(N,wc)和B=fir1(N,wc,win)表示設計一個具有線性相位的N階FIR低通或帶通濾波器，截止頻率wc為對奈奎斯特頻率進行歸一化的值。如果是帶通，wc為低端和高端截止頻率組成的向量。窗函數(shù)用變量win指明，默認為漢明窗(hamming)。

如果是高通或帶阻，用B=firl(N,wc.′highorstop′)格式,輸入?yún)?shù)加字符串‘high’或‘stop’。

(2)函數(shù)fir2：用頻率抽樣法計算FIR濾波系數(shù)，調用格式有以下兩種：

B=fir2(N,F,A)

B=fir2(N,F,A,win)

用頻率抽樣法設計一個具有線性相位的N階FIR濾波器，其頻率響應由向量F和A指定，默認為漢明窗(hamming)。如果改用其他窗函數(shù)，可用變量win指明。

(3)界面fdatool：利用MATLAB的濾波器設計與分析圖形界面(FilterDesignandAnalysisTool)可以快速地設計各類濾波器。使用方法是在MATLAB命令窗輸入fdatool命令即可打開設計界面，設置相關選項和輸入相關參數(shù)，即可得到濾波器參數(shù)。6.2.4

DSPLIB庫中的FIR濾波函數(shù)

DSPLIB庫有幾個現(xiàn)成的FIR濾波函數(shù),可以直接調用。(1)fir：直接型實FIR濾波，調用格式為

fir(DATA*x,DATA*h,DATA*r,DATA*dbuffer,ushortnh,ushortnx)

算法為變量說明：

*x為輸入數(shù)據(jù)實向量x的指針，x為DATA數(shù)據(jù)(Q15格式數(shù)表示為整型數(shù)據(jù))，nx為向量大小，x的存儲需滿足循環(huán)尋址要求；

*h為濾波系數(shù)h的指針，h=［b0,b1,b2,…］為DATA數(shù)據(jù)，nh為向量大小，h的存儲需滿足循環(huán)尋址要求；

*r為濾波輸出r的指針，r為DATA數(shù)據(jù)；

*dbuffer為濾波計算緩沖區(qū)，dbuffer的存儲需滿足循環(huán)尋址要求，第一次濾波時需初始化為0。

(2)firdec：直接型帶輸入抽取運算的實FIR濾波，調用格式為

firdec(DATA*x,DATA*h,DATA*r,DATA**dbuffer,ushortnh,ushortnx,ushortD)

算法為

D為抽取因子，其他變量說明均同fir函數(shù)。

(3)firinterp：直接型帶輸入插值運算的實FIR濾波，調用格式為

firinterp(DATA*x,DATA*h,DATA*r,DATA**dbuffer,ushortnh,ushortnx,ushortI)

算法為

I為抽取因子，其他變量說明均同fir函數(shù)。

(4)firs：系數(shù)對稱的實FIR濾波，調用格式為

firs(DATA*x,DATA*r,DATA**dbuffer,ushortnh2,ushortnx)

算法為

(5)firs2：通用系數(shù)對稱的實FIR濾波，調用格式為

firs2(DATA*x,DATA*h,DATA*r,DATA**dbuffer,ushortnh2,ushortnx)

算法為

同firs函數(shù)相比，系數(shù)向量不必存儲在程序存儲器中，適合多濾波器程序，但執(zhí)行效率較低。6.2.5調用DSPLIB庫文件的方法

設54xdsp.lib庫文件在C:＼ti＼c5400＼dsplib下，其相關頭文件在C:＼ti＼c5400＼dsplib＼include中,則CCS軟件中BuildOptions(構建選項)的設置方法如下。

(1)在Compiler選項卡的Preprocessor中的IncludeSearchPath(-i)中輸入頭文件地址,如圖6.8所示。圖6.8設置包含文件搜索途徑

(2)在Compiler選項卡的Basic中的LibarySearchPath(-i)中輸入庫文件地址，IncludeLibraries(-l)中輸入工程所包含的庫文件名,如圖6.9所示。

如果所建工程正在庫文件所在的文件夾，上述設置可以簡化。

在Preprocessor中的IncludeSearchPath(-i)中輸入引號中的符號“..＼..＼”即可，同樣，在Basic中的LibarySearchPath(-i)中也是輸入“..＼..＼”。圖6.9設置庫文件搜索途徑及包含的庫文件名6.3

IIR濾波器設計

6.3.1

IIR濾波器傳遞函數(shù)及實現(xiàn)結構

一個IIR數(shù)字濾波器的系統(tǒng)函數(shù)為(6.17)式中,一般有M≤N，稱為N階IIR系統(tǒng)。令z=ejω，可得其頻域特性為據(jù)此可求出其幅度特性和相位特性。(6.1８)式(6.18)的IIR濾波器時域輸入輸出關系為(6.19)

IR濾波器有多種實現(xiàn)結構，圖6.10(a)所示即為一個二階IIR濾波器直接I型結構，直觀明了，但需要4個存儲單元存儲x和y過去的值。若需要減少存儲單元，可使用直接II型結構，如圖6.10(b)所示。對于一個二階IIR濾波器來說，只需2個存儲單元。圖6.10二階IIR濾波器的直接I型和II型結構若要避免直接型結構的缺點，可采用級聯(lián)型或并聯(lián)型結構,如圖6.11所示。它將系統(tǒng)函數(shù)分解為若干個二階系統(tǒng)函數(shù)(稱為一節(jié))并以級聯(lián)或并聯(lián)的形式實現(xiàn)，這里把一階系統(tǒng)函數(shù)看做是二階的特例。由于各節(jié)結構相同，容易實現(xiàn)，同時又相互獨立，因此系數(shù)調整對其他各階沒有影響。圖6.11

IIR濾波器的結構6.3.2

IIR濾波系數(shù)的MATLAB輔助設計

例6.5設計一個數(shù)字切貝雪夫帶通濾波器，給定指標為(1)帶內波動δ1≤2dB，當200Hz≤f≤400Hz；

(2)帶外衰減δ2≥20dB，當f≤100Hz或f≥600Hz；

(3)抽樣頻率fs=2kHz。解：(1)求傳遞函數(shù)。

在MATLAB輸入：

fs=8e3;

%fs=8kHz

wp=500/(fs/2);

%Hz

ws=1000/(fs/2);%Hz

rp=3;

rs=50;

［n,wn］=buttord(wp,ws,rp,rs);

［b,a］=butter(n,wn);

freqz(b,a,512);畫出頻率特性運行MATLAB后，可以求出所設計的系統(tǒng)函數(shù)為(6.20)其中其頻率特性如圖6.12所示，其中橫坐標為歸一化頻率，乘以奈奎斯特頻率(fs/2)即為實際頻率，可以驗證其是否符合設計要求。圖6.12例6.5的IIR濾波器頻率特性

(2)系數(shù)歸一化。要解決DSP的實現(xiàn)問題最簡單的方法

就是直接采用浮點運算，但其運算量較大。C54x中，采

用Q15數(shù)據(jù)格式作乘加運算，可大大提高運算速度。為此，可考慮將輸入x(n)、輸出y(n)及系數(shù)向量b和a都統(tǒng)一為Q15

格式。因此，對式(6.20)的傳遞函數(shù)作適當變換，分子分母都除以略大于分子分母絕對值最大的數(shù),用MATLAB命令處理如下：

MAX=max(abs(［b,a］))*1.01；

b2=b/MAX;

a2=a/MAX;

處理后的系數(shù)向量b2和a2的每個系數(shù)絕對值都小于1，可用Q15表示。當輸入x(n)和系數(shù)向量[WTBZ]b2和a2[WTBZ]都用Q15格式數(shù)據(jù)表示后，濾波后的輸出y(n)是否能保證用Q15格式數(shù)據(jù)表示，還需要研究，可用MATLAB進行仿真驗證。有關MATLAB代碼如下：

fc=400;fs=fc*8

n=0:fs*8;

x=sin(2*pi*fc/fs*n)/2;

y2=filter(b2,a2,x);用調整后的系數(shù)向量b2和a2計算輸出。如果y(n)的每個絕對值都小于或等于1，則沒有溢出;如果有溢出，考慮到IIR系統(tǒng)是線性系統(tǒng)，可將輸入x(n)都統(tǒng)一縮小一定的幅度，以保證輸出不溢出，最簡單的方法就是將x(n)的值縮小一半，x(n)算術右移一位。

(3)轉化為DSP數(shù)據(jù)格式。由于CCS開發(fā)環(huán)境的匯編器

不接受小數(shù)數(shù)據(jù)，需將Q15格式數(shù)據(jù)轉化為整型數(shù)據(jù)。用MATLAB生成一個系數(shù)表，復制到一個空的.inc文件中。

B=round(b2*2∧15);A=round(a2*2∧15);

(4)必要時進行濾波器模型變換。使用級聯(lián)型IIR結構，可以方便、獨立地調整傳遞函數(shù)的零極點，且所用的存儲單元最少。在MATLAB中，將傳遞函數(shù)轉化為二階分式級聯(lián)模型的命令為

［sos,g］=tf2sos(num,den);式中，num、den分別為傳遞函數(shù)分子、分母多項式的系數(shù)，sos為二階分式級聯(lián)模型的系數(shù)，如下所示：對式(6.20)進行函數(shù)傳遞，可求出二階分式級聯(lián)模型的系數(shù)為同上,對sos進行歸一化，即可得到DSP編程所需的系數(shù)。6.3.3

DSPLIB庫的IIR濾波器程序

根據(jù)算法，可以自編匯編程序實現(xiàn)IIR濾波，也可以調用DSPLIB庫的IIR濾波函數(shù)來實現(xiàn)。IIR濾波函數(shù)主要有以下幾個函數(shù)。

(1)iircas4：二階級聯(lián)的IIR濾波函數(shù)，每節(jié)濾波器采用直接II型，4個系數(shù)。

shortiircas4(DATA*x,DATA*h,DATA*r,DATA**dbuffer,

ushortnbiq,ushortnx)

(2)iircas5：二階級聯(lián)的IIR濾波函數(shù)，每節(jié)濾波器采用直接II型，5個系數(shù)。

shortiircas5(DATA*x,DATA*h,DATA*r,DATA**dbuffer,ushortnbiq,ushortnx)

(3)iircas51：二階級聯(lián)的IIR濾波函數(shù)，每節(jié)濾波器采用直接I型，5個系數(shù)。

shortiircas51(DATA*x,DATA*h,DATA*r,DATA**dbuffer,ushortnbiq,ushortnx)

(4)iir32：二階級聯(lián)的IIR濾波函數(shù)，32位數(shù)據(jù)，每節(jié)濾波器采用直接II型，5個系數(shù)。

iir32(DATA*x,LDATA*h,DATA*r,LDATA**dbuffer,ushortnbiq,ushortnx)

6.4快速傅立葉變換的DSP實現(xiàn)

6.4.1基2復數(shù)FFT算法的原理

有限長序列的頻譜是周期性的連續(xù)函數(shù)，在數(shù)字頻率的一個周期ω=0~2π內，頻譜有無數(shù)個值。用DFT(離散傅立葉變換)可計算在一個數(shù)字頻率周期內的N個抽樣值，其計算公式為(6.21)式中,數(shù)字抽樣頻點為

FFT是離散傅立葉變換(DFT)的一種快速算法。由于我們在計算DFT時,一次復數(shù)乘法需用四次實數(shù)乘法和兩次實數(shù)加法,一次復數(shù)加法則需兩次實數(shù)加法。這樣，每計算一個X(k)需要4N次復數(shù)乘法和2N+2(N－1)=2(2N－1)次實數(shù)加法,因此整個DFT運算總共需要4N2次實數(shù)乘法和N×2(2N－1)=4N2－2N次實數(shù)加法。設復數(shù)序列x(n),n=0,1,…,N－1的長度為N=2L，L為整數(shù)，若長度N不是2L形式，則用數(shù)字0補足長度。將N=2L的復數(shù)序列x(n)按N的奇偶性分成x1(m)和x2(m)兩組m=0,1,…,

N/2－1，分別作DFT變換，得X1(k)和X2(k)，也就是說將一個N點的DFT分解成兩個N/2點的DFT，然后按式(6.22)重新組合成一個N點的DFT。(6.22)(6.2３)同直接計算DFT的運算量相比，將一個N點的DFT分解成兩個N/2點DFT，其運算量大約減少一半。如果繼續(xù)這個過程，運算量將以指數(shù)級下降。一個N點的DFT運算最終分解成N/2個2點的DFT運算，整個運算只需求（N／2）log2N次復數(shù)乘法運算，運算量比原來的N2次的復數(shù)乘法次數(shù)大大降低了(特別是N較大時)，這種算法稱為基2復數(shù)FFT算法。基2復數(shù)FFT算法將一個N點DFT運算每一級都分解成N/2個2點DFT運算，共有l(wèi)og2N級，但每級的加權系數(shù)WkN的個數(shù)不同。在log2N級(最后一級)，共有N/2個系數(shù)WkN,0≤k≤N/

2－1，每往前一級，系數(shù)個數(shù)減少一半，但系數(shù)的指數(shù)為前一級系數(shù)指數(shù)的2倍。第一級時，系數(shù)只有1個，為W0N=1，這時復數(shù)乘法就可以省去。6.4.2實序列FFT算法(RFFT)原理

當輸入為純實數(shù)的時候，可以利用左右對稱的特性，更好地計算DFT，稱這樣的RFFT優(yōu)化算法為包裝算法：首先將2N點實序列的連續(xù)輸入包裝成復數(shù)，稱為“進包”，然后運行N點復數(shù)FFT，最后將N點復數(shù)FFT輸出并合成2N點實序列的FFT輸出，稱為“開包裝”。

2N點實序列FFT算法如下：(6.24)把x(n)分成偶數(shù)點和奇數(shù)點，可得到(6.25)

n為偶數(shù)時，令n=2r;n為奇數(shù)時，令n=2r+1，則有(6.26)將X(k)分解為實部和虛部，即

X(k)=XR(k)+jXI(k)

(6.27)式中(6.28)(6.29)由FFT的性質可知道，實序列的DFT是共軛對稱的，即

X(k)=X*(2N－k)

(6.30)

或XR(k)=XR(2N－k),XI(k)=－XI(2N－k)。

下面通過N點復數(shù)FFT運算得到上面的X(k)。令N點復序列為y(r)，通過如下方法構造

y(r)=x(2r)+jx(2r+1)，r=0,1,2,…,N－1(6.31)則y(r)的DFTY(k)為(6.32)將Y(k)分解為實部和虛部，得