神經網絡第二章

1、會計學1神經網絡第二章神經網絡第二章2MATLAB是美國MathWorks公司開發(fā)的一種功能極其強大的高技術計算語言和內容極其豐富的軟件庫。它以矩陣和向量的運算以及運算結果的可視化為基礎，把廣泛應用于各個學科領域的數(shù)值分析、矩陣計算、函數(shù)生成、信號、圖形及圖象處理、建模與仿真等諸多強大功能集成在一個便于用戶使用的交互式環(huán)境之中，為使用者提供了一個高效的編程工具及豐富的算法資源。 關于MATLAB3 MATLAB與信號處理直接有關的工具箱與信號處理直接有關的工具箱 （ Toolbox） Signal Processing (信號處理工具箱） Wavelet (小波工具箱) Image Proce

2、ssing（圖象處理工具箱） Higher-Order Spectral Analysis （高階譜分析工具箱）4 與信號處理間接有關的工具箱：與信號處理間接有關的工具箱：Control System (控制系統(tǒng))Communication （通信）System Identification （系統(tǒng)辨識）Statistics （統(tǒng)計）Neural Network (神經網絡)5例例： z=peaks; surf(z)；61. rand.m 用來產生均值為0.5、幅度在 0 1 之 間 均 勻 分 布 的 偽 白 噪 聲 : u=rand(N,1) (rand(N)生成N階矩陣) 2122011

3、( ),12Nuuunu nN方差：如何改變 的方差( )u n( )u n與第二章內容有關的MATLAB文件方差函數(shù)var(u)標準差函數(shù)std(u)71. rand.m 用來產生均值為0.5、幅度在 0 1 之 間 均 勻 分 布 的 偽 白 噪 聲 : u=rand(N,1) (rand(N)生成N階矩陣) 2. randn.m 用來產生均值為零、方差為1 服從高斯（正態(tài)）分布的白噪聲信號 u=randn(1, N) ( )u n與第二章內容有關的MATLAB文件x=randn(1000,1)y=randn(1000,1)v=var(x)h=std(y)83.sinc ：用來產生 “si

4、nc” 函數(shù)：sinc函數(shù)定義為： 10sin ( )0sin()tc ttkttt為其它t=-4:0.1:4;x4=sinc(t); %產生抽樣函數(shù)plot(t,x4) 94. conv.m 用來實現(xiàn)兩個離散序列的線性卷積。其調用格式是：y=conv(x,h).若x(n)和y(n)的長度分別為M和N， 則返回值是長度為M+N-1的序列。例例 x(n)=3 4 5; h(n)=2 6 7 8，求其線性卷積。 MATLAB語句如下： x=3 4 5; h=2 6 7 8; y=conv(x,h)運行結果: y= 6 26 55 82 67 40 u兩序列的相關運算兩序列的相關運算12()( )(

5、)ny mxn xnmMATLAB實現(xiàn)：實現(xiàn)：y=xcorr(x1,x2)。x=3 4 5;h=2 6 7 8;y=xcorr(x,h)y =24 53 86 65 38 10 -010 5xcorr: 其互相關和自相關。格式是：(1)rxy=xcorr(x,y)：求x,y的互相關； (2)rx=xcorr(x,M,flag）：求x的自相關，M：rx的單邊長度，總長度為2M+1；flag是定標標志，若 flag=biased， 則表示是“有偏”估計，需將rx(m)都除以N，若flag=unbiased,則表示是“無偏”估計，需將rx(m)都除以（Nabs(m)）；若flag缺省，則rx不定標。

6、M和flag同樣適用于求互相關。1112u 例例 用用MATLAB求出離散序列求出離散序列 的的Z變換變換MATLAB程序如下：程序如下：syms k zf=0.5k; %定義離散信號定義離散信號Fz=ztrans(f) %對離散信號進行對離散信號進行Z變換變換u 運行結果如下：運行結果如下： Fz = 2*z/(2*z-1)13u 例例：求序列的：求序列的Z變換變換.syms n hn=sym(kroneckerDelta(n, 1) + kroneckerDelta(n, 2)+ kroneckerDelta(n, 3)Hz=ztrans(hn)Hz=simplify(Hz)u 運行結果如

7、下：運行結果如下：Fz= (z2 + z + 1)/z3( )(1)(4)f kkt2( )21zF zz14 1filter.m本文件用來求離散系統(tǒng)的輸出y(n) 。若系統(tǒng)的h(n)已知，由y(n)=x(n)*h(n)，用conv.m文件可求出y(n) 。 y=conv (x, h)filter文件是在A(z)、B(z)已知，但不知道h(n)的情況下求y(n)的。 調用格式是: y=filter(b, a, x) x, y, a 和 b都是向量。與逆與逆Z變換變換 相關的相關的matlab 函數(shù)函數(shù)15與逆與逆Z變換變換 相關的相關的matlab 函數(shù)函數(shù)2.impz.m在 A（z）、B（z

8、）已知情況下, 求系統(tǒng)的單位抽樣響應 h(n)。調用格式是： h = impz(b, a, N) 或 h,t=impz(b,a,N) N是所需的的長度。前者繪圖時n從1開始，而后者從0開始。16 3. residuez.m將X(z) 的有理分式分解成簡單有理分式的和，因此可用來求逆變換。調用格式： r,p,k= residuez(b,a)假如知道了向量r, p和k，利用residuez.m還可反過來求出多項式A(z)、B(z)。格式是 b,a= residuez(r,p,k)。 17頻率響應函數(shù)：頻率響應函數(shù)：freqz.m已知A(z)、B(z), 求系統(tǒng)的頻率響應?；镜恼{用格式是： H,w

9、=freqz(b,a,N,whole,Fs)N是頻率軸的分點數(shù)，建議N為2的整次冪；w是返回頻率軸座標向量，繪圖用；Fs是抽樣頻率，若Fs1，頻率軸給出歸一化頻率；whole指定計算的頻率范圍是從0FS,缺省時是從0FS/2.( )( )( )B zH zA z幅頻率響應：幅頻率響應：Hr=abs(H);相頻響應：相頻響應： Hphase=angle(H);解卷繞：解卷繞：Hphase=unwrap(Hphase);18下面幾個文件用于轉移函數(shù)與極零點之 間的相互轉換及極零點的排序：（1) tf2zpk.m, 調用z,p,k=tf2zpk(b,a) (2) zp2tf.m, 調用 b,a=zp

10、2tf (z,p,k) (3)roots.m， 調用 Z1=roots(b) (4) poly.m, 調用b =poly (Z1)與極零點有關的與極零點有關的MATLAB函數(shù)函數(shù)19 顯示離散系統(tǒng)的極零圖顯示離散系統(tǒng)的極零圖函數(shù)函數(shù)：zplane.m本文件可用來顯示離散系統(tǒng)的極零圖。其調用格式是： zplane(z,p), 或 zplane(b,a),前者是在已知系統(tǒng)零點的列向量列向量z和極點的列向量列向量p的情況下畫出極零圖，后者是在僅已知A(z)、B(z) 的情況下畫出極零圖。20-1-0.500.511.52-1-0.500.5122Real PartImaginary Part1221

11、222( )1 44(44)/(2) /H zzzzzzzz 用極零分析大致畫出幅頻響應及相頻響應：( )( )4 (1)4 (2)y nx nx nx n例1： 系統(tǒng)解：轉移函數(shù)：b=1 -4 4;a=1;z,p,k=tf2zpk(b,a)zplane(b,a)zplane(z,p)r,p,k= residuez(b,a)b,a= residuez(r,p,k)z=2; 2p=0; 0K=1r =; p=;k=1 -4 4;211221222( )1 44(44)/(2) /H zzzzzzzz 用極零分析大致畫出幅頻響應及相頻響應：( )( )4 (1)4 (2)y nx nx nx n例

12、1： 系統(tǒng)解：轉移函數(shù)：b=1 -4 4;a=1;H,w=freqz(b,a,128,whole,1)Hr=abs(H);Hphase=angle(H);Hphaseu=unwrap(Hphase);subplot(311),plot(Hr)subplot(312),plot(Hphase)subplot(313),plot(Hphaseu)00.20.40.60.810510 Amplitude Freq. Res.00.20.40.60.81-505 wrap Phase Res.00.20.40.60.81-15-10-50unwrap Phase Res.22|()| 1jH e11(

13、 )H zzz例例2：( ),02 相位的卷繞相位的卷繞 (wrapping) 解卷繞解卷繞 b=1;a=0,1;H,w=freqz(b,a,256,whole,1);Hr=abs(H);Hphase=angle(H);Hphase1=unwrap(Hphase); 050100150200250300111111050100150200250300-4-20240501001502002503000246823-1-2-3-4-1-2-3-4.1836+.7344z +1.1016z +.7374z +.183116z( )1-3.0544z +3.8291z -2.2925z +.5507

14、5z00H z 例：例： 給定系統(tǒng)給定系統(tǒng)求：極零圖 單位抽樣響應 頻率響應H,w=freqz(b,a,256,whole,1); Hr=abs(H); Hphase=angle(H); Hphase=unwrap(Hphase); h,t=impz(b,a,40);stem(t,h,.);grid on;zplane(b,a);b=.1836 .7344 1.1016 .7374 .1836/100a =1 -3.0544 3.8291 -2.2925 .55075解：24-1.5-1-0.500.51-1-0.8-0.6-0.4-0.200.20.40.60.81Real PartImag

15、inary Part極零圖極零圖 zplane(b,a); 250510152025303540-0.1-0.0500.050.10.150.20.25單位抽樣響應h,t=impz(b,a,40);stem(t,h,.);grid on;26頻率響應00.10.20.30.40.50.60.70.80.9100.511.5 Amplitude Freq. Res.00.10.20.30.40.50.60.70.80.91-505 wrap Phase Res.00.10.20.30.40.50.60.70.80.91-20-100unwrap Phase Res.Hphase=angle(H)

16、; Hphaseu=unwrap(Hphase); H,w=freqz(b,a,256,whole,1); Hr=abs(H); subplot(311),plot(Hr)subplot(312),plot(Hphase)subplot(313),plot(Hphaseu)27下面幾個文件用于轉移函數(shù)、極零點與二階子系統(tǒng)二階子系統(tǒng)sos(Second-Order Section)級聯(lián)之間的相互轉換：（1) tf2sos.m, 調用sos,G=tf2sos(b,a) (2) sos2tf.m, 調用 b,a=sos2tf (sos,G) (3) sos2zp.m, 調用z,p,k= sos2zp

17、 (sos,G) (4) zp2sos.m, 調用 sos,G=zp2sos (z,p,k)與信號流圖有關的與信號流圖有關的MATLAB函數(shù)函數(shù)12,0,1,212,1,2( ),1,1kkkkkkzzHzkLazaz,0,1,2,1,21,1,kkkkkaakLsos是一是一L6的矩陣，每行元素如下排列：的矩陣，每行元素如下排列：281buttord.m 確定 LP DF、或 LP AF的階次；（1） N, Wn = buttord(Wp, Ws, Rp, Rs)；（2）N, Wn = buttord(Wp, Ws, Rp, Rs,s)：與本章內容有關的MATLAB文件(1)對應 數(shù)字濾波器

18、。其中 Wp, Ws分別是通帶和阻帶的截止頻率，其值在 01 之間，1對應抽樣頻率的一半(歸一化頻率)。對低通和高通，Wp, Ws都是標量，對帶通和帶阻，Wp, Ws是12的向量。Rp, Rs 分別是通帶和阻帶的衰減(dB)。N是求出的相應低通濾波器的階次，Wn是求出的3dB頻率，它和Wp稍有不同。(2)對應模擬濾波器，各變量含意和（1）相同，但Wp, Ws及Wn的單位為弧度/秒，它們實際上是頻率。292buttap.m 設計模擬低通(Butt)原型濾波器。 z, p, k=buttap(N): N是欲設計的低通原型濾波器的階次，z, p, k是設計出的極點、零點及增益。3lp2lp.m、lp

19、2hp.m、lp2bp.m, lp2bs.m將模擬低通原型轉換為實際的低通、高通、帶通及帶阻濾波器。b, a 是AF LP 的分子、分母的系數(shù)向量，B, A是轉換后的的分子、分母的系數(shù)向量；(1)中，Wo是低通或高通濾波器的截止頻率；B, A=lp2lp(b, a, Wo) B, A=lp2hp(b, a, Wo)(1)B, A=lp2bp(b, a, Wo , Bw) B, A=lp2bs(b, a, Wo , Bw)(2)(2)中Wo是帶通或帶阻濾波器中心頻率，Bw是其帶寬。304bilinear.m ：雙線性變換，由模擬濾波器 得到數(shù)字濾波器。 Bz, Az=bilinear(B, A,

20、 Fs) 式中B, A分別是G（s）的分子、分母多項式的系數(shù)向量，Bz, Az分別是H(z)的分子、分母多項式的系數(shù)向量，F(xiàn)s是抽樣頻率。315butter.m 用來直接設計Butterworth數(shù)字濾波器，實際上它把 buttord.m, buttap.m, lp2lp.m, bilinear.m等文件都包含了進去，從而使設計過程更簡捷。格式(1)(3) 用來設計數(shù)字濾波器，B，A分別是H(z)的分子、分母多項式的系數(shù)向量，Wn是通帶截止頻率，范圍在01之間。若Wn是標量，(1)用來設計低通數(shù)字濾波器，若Wn是12的向量，則(1) 用來設計數(shù)字帶通濾波器； (2)用來設計數(shù)字高通濾波器； (

21、3) 用來設計數(shù)字帶阻濾波器，顯然，這時的Wn是12的向量；格式(4) 用來設計模擬濾波器。(1) B,A=butter(N,Wn); (2) B,A=butter(N,Wn,high); (3) B,A=butter(N,Wn,stop); (4) B,A=butter(N,Wn,s)32nh,w=freqz(bz,az,256,Fs*1000);nplot(w,abs(h);grid on;設計 IIR LP DF，3dB,20dBpssDF:100Hz,300Hz,1000HzpsffF33nh2,w2=freqz(bz,az,500,Fs);nplot(w1/2/pi,abs(h1),

22、w2,abs(h2),k);grid on;設計 IIR LP DF，3dB,20dBpssDF:100Hz,300Hz,1000HzpsffF34設計 IIR LP DF，3dB,20dBpssDF:100Hz,300Hz,1000HzpsffF35非線性關系DF:100Hz,300Hz,1000Hz0.2 ,0.6 ,2AF:tan(/2)685.82109(Hz)2tan(/2)2452.762438(Hz)psspsppssssffFTT 設計的 AF 并不是按給定的技術指標，但當再由 變回 后，保證了 DF的技術要求。sz又稱為頻率的預變形（Freq. Warping) 。例如 ：

23、2tan(/2)sT 抽樣頻率pf .2sf .236 給出給出數(shù)字高通數(shù)字高通的技術要求的技術要求)(zHdhpspsp, 得到得到模擬高通模擬高通的技術要求的技術要求)(sHahpspsp,)2tan(1step得到得到模擬低通模擬低通的技術要求的技術要求)(pGspsp,p12step設計出設計出)(pG3stepspp4step 得到得到模擬高通轉移模擬高通轉移函數(shù)函數(shù))(sHahp11zzs5step最后得到最后得到數(shù)字高通轉移數(shù)字高通轉移函數(shù)函數(shù))(zHdhp數(shù)字高通濾波器設計步驟7.6 數(shù)字高通, 帶通及帶阻濾波器的設計psqpq/137對 帶通（BP）、帶阻（BS）數(shù)字濾波器的

24、設計，只需改變圖中 Step2 和 Step4：22231BWBW 帶阻31213()sps 22231BWBW 帶通21331()sps 382N 121)(2pppG) 1() 1() 1(22222zzzpBW432142641. 0307. 1237. 2637. 11)21 (0201. 0)(zzzzzzzHdbps3dB,18dB,F2000Hzps要求：按上述轉換辦法，可以求出：例6.6.2： 設計一 IIR BP DF,要求：通帶頻率范圍 ： 300Hz 400Hz ;阻帶頻率范圍 ：200Hz、500Hz31 ffslshff391buttord.m 確定 LP DF、或

25、LP AF的階次；（1） N, Wn = buttord(Wp, Ws, Rp, Rs)；（2）N, Wn = buttord(Wp, Ws, Rp, Rs,s)：與本章內容有關的MATLAB文件(1)對應 數(shù)字濾波器。其中 Wp, Ws分別是通帶和阻帶的截止頻率，其值在 01 之間，1對應抽樣頻率的一半(歸一化頻率)。對低通和高通，Wp, Ws都是標量，對帶通和帶阻，Wp, Ws是12的向量。Rp, Rs 分別是通帶和阻帶的衰減(dB)。N是求出的相應低通濾波器的階次，Wn是求出的3dB頻率，它和Wp稍有不同。(2)對應模擬濾波器，各變量含意和（1）相同，但Wp, Ws及Wn的單位為弧度/秒

26、，它們實際上是頻率。402buttap.m 設計模擬低通(Butt)原型濾波器。 z, p, k=buttap(N): N是欲設計的低通原型濾波器的階次，z, p, k是設計出的極點、零點及增益。3lp2lp.m、lp2hp.m、lp2bp.m, lp2bs.m將模擬低通原型轉換為實際的低通、高通、帶通及帶阻濾波器。b, a 是AF LP 的分子、分母的系數(shù)向量，B, A是轉換后的的分子、分母的系數(shù)向量；(1)中，Wo是低通或高通濾波器的截止頻率；B, A=lp2lp(b, a, Wo) B, A=lp2hp(b, a, Wo)(1)B, A=lp2bp(b, a, Wo , Bw) B, A

27、=lp2bs(b, a, Wo , Bw)(2)(2)中Wo是帶通或帶阻濾波器中心頻率，Bw是其帶寬。414bilinear.m ：雙線性變換，由模擬濾波器 得到數(shù)字濾波器。 Bz, Az=bilinear(B, A, Fs) 式中B, A分別是G（s）的分子、分母多項式的系數(shù)向量，Bz, Az分別是H(z)的分子、分母多項式的系數(shù)向量，F(xiàn)s是抽樣頻率。425butter.m 用來直接設計Butterworth數(shù)字濾波器，實際上它把 buttord.m, buttap.m, lp2lp.m, bilinear.m等文件都包含了進去，從而使設計過程更簡捷。格式(1)(3) 用來設計數(shù)字濾波器，B

28、，A分別是H(z)的分子、分母多項式的系數(shù)向量，Wn是通帶截止頻率，范圍在01之間。若Wn是標量，(1)用來設計低通數(shù)字濾波器，若Wn是12的向量，則(1) 用來設計數(shù)字帶通濾波器； (2)用來設計數(shù)字高通濾波器； (3) 用來設計數(shù)字帶阻濾波器，顯然，這時的Wn是12的向量；格式(4) 用來設計模擬濾波器。(1) B,A=butter(N,Wn); (2) B,A=butter(N,Wn,high); (3) B,A=butter(N,Wn,stop); (4) B,A=butter(N,Wn,s)436cheb1ord.m 求Cheb-型濾波器的階；7cheb1ap.m 設計原型低Cheb

29、-I型模擬濾波器；8cheby1.m 直接設計數(shù)字Cheb-濾波器。以上三個文件的調用格式和對應的Butterworth濾波器的文件類似。449cheb2ord.m; 10. ellipord.m; 11.cheb2ap.m; 12. ellipap.m; 13.besselap.m; 14. cheby2.m; 15. ellip.m; 16.besself.m17impinvar.m 用沖激響應不變法實現(xiàn)頻率轉換； 對應 Cheby-II、橢圓 IIR 濾波器45產生窗函數(shù)的文件有八個：1. bartlett（三角窗）; 2. blackman（布萊克曼窗） ; 3. boxcar（矩形窗

30、）; 4. hamming（哈明窗）; 5. hanning（漢寧窗）; 6. triang（三角窗）;7. chebwin（切比雪夫窗）; 8 .kaiser（凱賽窗）; 兩端為零兩端不為零調用方式都非常簡單請見help文件稍為復雜469fir1.m 用“窗函數(shù)法”設計FIR DF。調用格式： （1）b = fir1(N,Wn); (2) b = fir1(N,Wn,high); (3) b = fir1(N,Wn, stop); N：階次，濾波器長度為N1；Wn：通帶截止頻率，其值在01之間，1對應 Fs/2； b: 濾波器系數(shù)。 格式（2）用來設計高通濾波器， 格式（3）用來設計帶阻濾波

31、器。 格式（1），若Wn為標量，則設計低通濾波器，若Wn是12的向量，則用來設計帶通濾波器，若Wn是1L的向量，則可用來設計L帶濾波器。47對格式（1），若Wn為標量，則設計低通濾波器，若Wn是12的向量，則用來設計帶通濾波器，若Wn是1L的向量，則可用來設計L帶濾波器。這時，格式（1）要改為: b = fir1(N,Wn, DC-1), 或 b = fir1(N,Wn, DC-0)前者保證第一個帶為通帶，后者保證第一個帶為阻帶。在上述所有格式中，若不指定窗函數(shù)的類型，fir1自動選擇Hamming窗。指定窗函數(shù)格式： （4）b = fir1(N,Wn,wind); 例 b = fir1(N,

32、Wn,boxcar(N+1); 指定矩形窗4810fir2.m 本文件采用“窗函數(shù)法”設計具有任意幅 頻相應的FIR 數(shù)字濾波器。其調用格式是： b = fir2(N, F, M); F是頻率向量，其值在01之間，M是和F相對應 的所希望的幅頻相應。如同fir1, 缺省時自動選用 Hamming窗。例 ：設計一多帶濾波器，要求頻率在0.20.3, 0.60.8 之間為1，其余處為零。 設計結果如下：4905101520253035-0.500.5020406080100-0.500.500.10.20.30.40.500.51N=30,90時幅頻響應響應及理想幅頻響應；N=30N=90( )h

33、 n()jH e5013 firls.m 用最小平方法設計線性相位FIR濾波器，可設計任意給定的理想幅頻響應；14 fircls.m用帶約束的最小平方法設計線性相位FIR濾波器，可設計任意給定的理想幅頻響應；15 fircls1.m 用帶約束的最小平方方法設計線性相位FIR低通和高通濾波器。16 sgolay.m 用來設計 Savitzky-Golay FIR 平滑濾波器，其原理見9.1.1節(jié) 17 firrcos.m 用來設計低通線性相位FIR濾波器，其過渡帶為余弦函數(shù)形狀。519fir1.m 用“窗函數(shù)法”設計FIR DF。調用格式： （1）b = fir1(N,Wn); (2) b =

34、fir1(N,Wn,high); (3) b = fir1(N,Wn, stop); N：階次，濾波器長度為N1；Wn：通帶截止頻率，其值在01之間，1對應 Fs/2； b: 濾波器系數(shù)。 格式（2）用來設計高通濾波器， 格式（3）用來設計帶阻濾波器。 格式（1），若Wn為標量，則設計低通濾波器，若Wn是12的向量，則用來設計帶通濾波器，若Wn是1L的向量，則可用來設計L帶濾波器。52對格式（1），若Wn為標量，則設計低通濾波器，若Wn是12的向量，則用來設計帶通濾波器，若Wn是1L的向量，則可用來設計L帶濾波器。這時，格式（1）要改為: b = fir1(N,Wn, DC-1), 或 b =

35、 fir1(N,Wn, DC-0)前者保證第一個帶為通帶，后者保證第一個帶為阻帶。在上述所有格式中，若不指定窗函數(shù)的類型，fir1自動選擇Hamming窗。指定窗函數(shù)格式： （4）b = fir1(N,Wn,wind); 例 b = fir1(N,Wn,boxcar(N+1); 指定矩形窗53例7.1.1.設計低通 DF FIR, 令截止頻率 0. 25, 取 M10, 20，40，觀察其幅頻響應的特點.cclear all;N=10;b1=fir1(N,0.25,boxcar(N+1); b3=fir1(2*N,0.25,boxcar(2*N+1); b4=fir1(4*N,0.25,box

36、car(4*N+1); M=128;h1=freqz(b1,1,M);h3=freqz(b3,1,M);h4=freqz(b4,1,M);f=0:0.5/M:0.5-0.5/M;plot(f,abs(h1),f,abs(h3),f,abs(h4);grid; axis(0 0.5 0 1.2)()jdHe22cc154例7.1.2： 理想差分器及其設計clear all;N=40;n=0:N;b1=fir1(N,0.25,boxcar(N+1); b2=fir1(N,0.25,hamming(N+1); win=hamming(N+1);for n=1:N+1 if (n-1-N/2)=0;

37、b1(n)=0; else b1(n)=(-1)(n-1-N/2)/(n-1-N/2); end endfor n=1:N+1 if (n-1-N/2)=0; b2(n)=0; elseb2(n)=win(n)*(-1)(n-1-N/2)/(n-1-N/2); end endM=128;h1=freqz(b1,1,M);h2=freqz(b2,1,M);% h=freqz(b,1,M);f=0:0.5/M:0.5-0.5/M;hd=2*pi*f;plot(f,abs(h1),f,abs(h2),f,hd,k-)/2( )(/ 2) ( )( 1)( )/ 2dn Mh nh nMw nw nn

38、M5511. remez.m 設計Chebyshev最佳一致逼近FIR濾波器、Hilbert變換器和差分器。調用格式是： (1) b=remez(N, F, A); (2) b=remez(N, F, A, W); （3）b=remez(N,F,A,W,Hilbert); (4) b=remez(N, F, A,W, differentiator)N是給定的濾波器的階次，b是設計的濾波器的系數(shù)，其長度為N1；F是頻率向量，A是對應F的各頻段上的理想幅頻響應，W是各頻段上的加權向量。56F、A及W的指定方式和例7.4.1和7.4.2所討論過的一樣，唯一的差別是F的范圍為01，而非00.5, 1對

39、應抽樣頻率的一半。需要指出的是，若b的長度為偶數(shù)，設計高通和帶阻濾波器時有可能出現(xiàn)錯誤，因此，最好保證b的長度為奇數(shù)，也即N應為偶數(shù)。57例7.4.1： 設計低通 FIR DF:0.6p,0.7sb=remez(N, F, A, W)0.60.71clear all;f=0 .6 .7 1;% 給定頻率軸分點；A=1 1 0 0;% 頻率分點上理想幅頻響應；weigh=1 10;% 頻率分點上的加權；b=remez(32,f,A,weigh);% 設計出切比雪夫最佳一致逼近濾波器;h,w=freqz(b,1,256,1);h=abs(h);h=20*log10(h);figure(1);ste

40、m(b,.);grid;figure(2);plot(w,h);grid; 調整通帶、阻帶的加權及濾波器的長度。調整N或W的結果58例7.4.2： 設計多阻帶濾波器，抽樣頻率500Hz, 在50Hz、 100Hz 及150Hz處陷波。 通帶加權為8，阻帶為1-17dB通帶、阻帶加權都是1-25dB59例7.4.2： 設計多阻帶濾波器，抽樣頻率500Hz, 在50Hz、 100Hz 及150Hz處陷波。 clear all;f=0 .14 .18 .22 .26 .34 .38 .42 .46 .54 .58 .62 .66 1;A=1 1 0 0 1 1 0 0 1 1 0 0 1 1;wei

41、gh1=8 1 8 1 8 1 8;b1=remez(64,f,A,weigh1); h1,w1=freqz(b1,1,256,1);h1=abs(h1);h1=20*log10(h1);subplot(211);plot(w1,h1);grid;axis(0 0.5 -60 10)title(N=32,weight=8 1 8 1 8 1 8,FontSize,14,Color,r)10.140.180.260.220.340.380.420.460.580.620.540.66250Hz6012remezord.m 本文件用來確定在用Chebyshev最佳一致逼近設計FIR濾波器時所需要的

42、濾波器階次。其調用格式是： N, Fo, Ao, W = remezord(F, A, DEV, Fs)。F、A的含意同文件remez，DEV是通帶和阻帶上的偏差；輸出的是適合要求的濾波器階次N、頻率向量Fo、幅度向量Ao和加權向量W。若設計者事先不能確定要設計的濾波器的階次，那么，調用remezord后，就可利用這一族參數(shù)調用remez, 即 b=remez(N, Fo, Ao, W)，從而設計出所需要濾波器。因此，remez和remezord常結合起來使用。需要說明的是，remezord給出的階次N有可能偏低，這時適當增加N即可；另外，最好判斷一下，若N為奇數(shù)，就令其加一，使其變?yōu)榕紨?shù)，這

43、樣b的長度為奇數(shù)。 61fftfilt.m 用疊接相加法實現(xiàn)長序列卷積。格式是: y=fftfilt(h,x) 或 y=fftfilt(h, x，N)與本章有關的 MATLAB 文件記 的長度為 ， 的長度為 。 若采用第一個調用方式，程序自動地確定對 分段的長度 及做FFT的長度 ， 顯然， 是最接近 的2的整次冪。分的段數(shù)為 。( )x nxN( )h nM( )x n()LMNL/xNLN采用第二個調用方式，使用者可自己指定做FFT的長度。建議使用第一個調用方式。62clear;% 用疊接相加法，計算濾波器系數(shù)用疊接相加法，計算濾波器系數(shù)h和輸入信號和輸入信號x的卷積的卷積% 其中其中h

44、為為10階階hanning窗，窗，x是帶有高斯白噪的正弦信號是帶有高斯白噪的正弦信號h=fir1(10,0.3,hanning(11);% h: is the impulse responseN=500;p=0.05;f=1/16; % of a low-pass filter.u=randn(1,N)*sqrt(p); % u:white noises=sin(2*pi*f*0:N-1); % s:sine signalx=u(1:N)+s; % x: a long sequence;y=fftfilt(h,x); % y=x*hsubplot(211)plot(x);subplot(212

45、)plot(y);例3.9.1令x(n)為一正弦加白噪聲信號，長度為500， h(n)是用fir1.m文件設計出的一個低通FIR濾波器，長度為11.試用fftfilt實現(xiàn)長序列的卷積實現(xiàn)長序列的卷積63050100150200250300350400450500-2-1012050100150200250300350400450500-2-101264clear;% 生成濾波器系數(shù)h和混有高斯白噪的正弦信號xh=fir1(10,0.3,hanning(11);N=500;p=0.05;f=1/16;u=randn(1,N)*sqrt(p);%s=sin(2*pi*f*0:N-1);x=u(1:

46、N)+s;% 將x分為長度為L的小段L=20;M=length(h);y=zeros(1,N+M-1);tempy=zeros(1,M+L-1);tempX=zeros(1,L);for k=0:N/L-1 tempx(1:L)=x(k*L+1:(k+1)*L); tempy=conv(tempx,h); y=y+zeros(1,k*L),tempy,zeros(1,N-(k+1)*L);endsubplot(211);plot(x)subplot(212);plot(y(1:N)65050100150200250300350400450500-2-1012050100150200250300

47、350400450500-2-101266 hilbert.m 文件用來計算信號Hilbert變換。調用的格式是: y=hilbert(x),y的實部就是 ，虛部是的Hilbert變換 。所以，y 實際上是 x 的解析信號。( )x n( )x n67 czt.m 調用格式是： Xczt(x, M, W, A) 。x是待變換的時域信號，其長度設為N，M是變換的長度，W確定變換的步長，A確定變換的起點。若M=N, A=1, 則CZT變成DFT。00,jjAeWeA=exp(j*2*pi*f0/fs);W=exp(-j*2*pi*DELf/fs);68例例。clear all; % 產生兩個正弦加

48、白噪聲； N=256; f1=.1;f2=.2;fs=1; a1=5;a2=3; w=2*pi/fs; x=a1*sin(w*f1*(0:N-1)+a2*sin(w*f2*(0:N-1)+randn(1,N); % 應用FFT 求頻譜； subplot(3,1,1); plot(x(1:N/4); f=-0.5:1/N:0.5-1/N; X=fft(x); y=ifft(X); subplot(3,1,2); plot(f,fftshift(abs(X);subplot(3,1,3); plot(real(y(1:N/4); 69010203040506070-10010-0.5-0.4-0.

49、3-0.2-0.100.10.20.30.40.50200400600010203040506070-1001070024681012141618200501000246810121416182005010077.588.599.51010.505010071 程序clear all；% 構造三個不同頻率的正弦信號的疊加作為試驗信號N=128;f1=8;f2=8.22;f3=9;fs=40;stepf=fs/N;n=0:N-1;t=2*pi*n/fs;n1=0:stepf:fs/2-stepf;x=sin(f1*t)+sin(f2*t)+sin(f3*t);M=N;W=exp(-j*2*pi/

50、M);% A=1時的czt變換A=1;Y1=czt(x,M,W,A);subplot(311)plot(n1,abs(Y1(1:N/2);grid on;72% DTFTY2=abs(fft(x);subplot(312)plot(n1,abs(Y2(1:N/2);grid on;% 詳細構造A后的czt M=60;f0=7.2;DELf=0.05;A=exp(j*2*pi*f0/fs);W=exp(-j*2*pi*DELf/fs);Y3=czt(x,M,W,A);n2=f0:DELf:f0+(M-1)*DELf;subplot(313);plot(n2,abs(Y3);grid on;731

51、filtfilt.m 本文件實現(xiàn)零相位濾波。其調用格式是：y=filtfilt(B, A, x) 。式中B是 的分子多項式，A是分母多項式，x是待濾波信號，y是濾波后的信號。 clear;N=32; n=-N/2:N+N/2; w=0.1*pi;x=cos(w*n)+cos(2*w*n);subplot(311);stem(n,x,.);grid on; xlabel(n);b=0.06745 0.1349 0.06745; a=1 -1.143 0.4128;y=filtfilt(b,a,x); % 用給定系統(tǒng)(b,a)對信號 x 作零相位濾波；y1=filter(b,a,x); % 用給定

52、系統(tǒng)(b,a)對信號 x 作低通濾波；subplot(312);stem(n,y,.);grid on; xlabel(n);subplot(313);stem(n,y1,.);grid on; xlabel(n);( )H z與本章內容有關的MATLAB文件74 2grpdelay.m 求系統(tǒng)的群延遲。調用格式 gd w=grpdelay(B, A, N) , 或 gd F=grpdelay(B, A, N, FS)式中B和A仍是 的分子、分母多項式，gd是群延遲，w、F是頻率分點，二者的維數(shù)均為N；FS為抽樣頻率，單位為Hz。( )H z00.511.522.533.500.511.522

53、.53 Amplitude Freq. Res.753deconv.m ：實現(xiàn)系統(tǒng)的反卷積，其調用格式： q,r=deconv(y,x); 也用來實現(xiàn)多項式除法。clear all;k=0:1:7; x=k+1; h=ones(1,4); y=conv(h,x); % y = x * h;q,r=deconv(y,x); % 由 y,x 作反卷積，求出 h；q1,r1=deconv(y,h); % 由 y,h 作反卷積，求出 x； subplot(321);stem(h,.b);ylabel( h(n);subplot(322);stem(x,.b);ylabel( x(n);subplot(

54、323);stem(y,.b);ylabel( y(n);subplot(325);stem(q,.r);ylabel( q(n);subplot(326);stem(q1,.r);ylabel( q1(n);clear all;% 實現(xiàn)多項式除法q=(x3+1)/(x+1)y=1 0 0 1; x=1 1; q,r=deconv(y,x); q763tf2latc.m 和latc2tf.m：實現(xiàn)轉移函數(shù)和Lattice 系數(shù)之間的相互轉換。tf2latc的調用格式是：(1) k=tf2latc(b), (2) k=tf2latc(1,a), (3) k, c=tf2latc(b,a), 其中（1）對應全零系統(tǒng)，（2）對應全極系統(tǒng)，（3）對應極零系統(tǒng)。latc2tf的調用格式和tf2latc正好相反。需要說明的是，tf2latc求出的Lattice系數(shù)k和本書求出的k差一個負號，這是由于我們在圖中用的是k。774. latcfilt.m 用來實現(xiàn)Lattice 結構下的信號濾波。調用格式是： (1) y, g=latcfilt(k, x)： 對應全零系統(tǒng) (2) y, g=latcfilt(k, 1, x)：對應全極系統(tǒng) (3) y, g=latcfilt(k, c, x)：對應極零系統(tǒng) x是待濾波的信號，y是用Lattice 結構作正向濾波的輸出，g是作反