神經(jīng)網(wǎng)絡(luò) 第二章 - 山東大學(xué)

1、2022-5-91Zhongguo Liu_Biomedical Engineering_Shandong Univ.Biomedical Signal processingmatlab 信號處理函數(shù)信號處理函數(shù)Zhongguo LiuBiomedical EngineeringSchool of Control Science and Engineering, Shandong UniversityMATLAB是美國MathWorks公司開發(fā)的一種功能極其強(qiáng)大的高技術(shù)計(jì)算語言和內(nèi)容極其豐富的軟件庫。它以矩陣和向量的運(yùn)算以及運(yùn)算結(jié)果的可視化為基礎(chǔ)，把廣泛應(yīng)用于各個(gè)學(xué)科領(lǐng)域的數(shù)值分析、矩陣計(jì)算、

2、函數(shù)生成、信號、圖形及圖象處理、建模與仿真等諸多強(qiáng)大功能集成在一個(gè)便于用戶使用的交互式環(huán)境之中，為使用者提供了一個(gè)高效的編程工具及豐富的算法資源。 關(guān)于MATLAB MATLAB與信號處理直接有關(guān)的工具箱與信號處理直接有關(guān)的工具箱 （ Toolbox） Signal Processing (信號處理工具箱） Wavelet (小波工具箱) Image Processing（圖象處理工具箱） Higher-Order Spectral Analysis （高階譜分析工具箱） 與信號處理間接有關(guān)的工具箱：與信號處理間接有關(guān)的工具箱：Control System (控制系統(tǒng))Communicatio

3、n （通信）System Identification （系統(tǒng)辨識）Statistics （統(tǒng)計(jì)）Neural Network (神經(jīng)網(wǎng)絡(luò))例例： z=peaks; surf(z)；1. rand.m 用來產(chǎn)生均值為0.5、幅度在 01之間均勻分布的偽白噪聲: u=rand(N,1) (rand(N)生成N階矩陣) 2122011( ),12Nuuunu nN方差：如何改變 的方差( )u n( )u n與第二章內(nèi)容有關(guān)的MATLAB文件方差函數(shù)var(u)標(biāo)準(zhǔn)差函數(shù)std(u)1. rand.m 用來產(chǎn)生均值為0.5、幅度在 01之間均勻分布的偽白噪聲: u=rand(N,1) (rand(

4、N)生成N階矩陣) randn.m 用來產(chǎn)生均值為零、方差為1 服從高斯（正態(tài)）分布的白噪聲信號 u=randn(1, N) ( )u n與第二章內(nèi)容有關(guān)的MATLAB文件x=randn(1000,1)y=randn(1000,1)v=var(x)h=std(y)3.sinc ：用來產(chǎn)生 “sinc” 函數(shù)：sinc函數(shù)定義為： 10sin ( )0sin()tc ttkttt為其它t=-4:0.1:4;x4=sinc(t); %產(chǎn)生抽樣函數(shù)plot(t,x4) 4. conv.m 用來實(shí)現(xiàn)兩個(gè)離散序列的線性卷積。其調(diào)用格式是：y=conv(x,h).若x(n)和y(n)的長度分別為M和N，

5、則返回值是長度為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ùn)行結(jié)果: y= 6 26 55 82 67 40 u兩序列的相關(guān)運(yùn)算兩序列的相關(guān)運(yùn)算12()( )()ny mxn xnmMATLAB實(shí)現(xiàn)：實(shí)現(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 -0 5xcorr: 其互相關(guān)和自相關(guān)。格式是：(1)rxy=xcorr(x,y)：求x,y的互相關(guān)； (2)rx=xcorr(

6、x,M,flag）：求x的自相關(guān)，M：rx的單邊長度，總長度為2M+1；flag是定標(biāo)標(biāo)志，若 flag=biased， 則表示是“有偏”估計(jì)，需將rx(m)都除以N，若flag=unbiased,則表示是“無偏”估計(jì)，需將rx(m)都除以（Nabs(m)）；若flag缺省，則rx不定標(biāo)。M和flag同樣適用于求互相關(guān)。第三章 Z變換. 在在MATLAB語言中有專門對信號進(jìn)行正反語言中有專門對信號進(jìn)行正反Z變換的函數(shù)變換的函數(shù)ztrans( ) 和和itrans( )。其調(diào)用格式分別如下：。其調(diào)用格式分別如下：uF=ztrans( f ) 對對f(n)進(jìn)行進(jìn)行Z變換，其結(jié)果為變換，其結(jié)果為F(

7、z)uF=ztrans(f,v) 對對f(n)進(jìn)行進(jìn)行Z變換，其結(jié)果為變換，其結(jié)果為F(v)uF=ztrans(f,u,v) 對對f(u)進(jìn)行進(jìn)行Z變換變換,其結(jié)果為其結(jié)果為F(v)uf=itrans ( F ) 對對F(z)進(jìn)行進(jìn)行Z反變換反變換,其結(jié)果為其結(jié)果為f(n)uf=itrans(F,u) 對對F(z)進(jìn)行進(jìn)行Z反變換，其結(jié)果為反變換，其結(jié)果為f(u)uf=itrans(F,v,u ) 對對F(v)進(jìn)行進(jìn)行Z反變換反變換,其結(jié)果為其結(jié)果為f(u)u注意：注意： 在調(diào)用函數(shù)在調(diào)用函數(shù)ztran( )及及iztran( )之前，要用之前，要用syms命令對所有需要用到的變量（如命令對所

8、有需要用到的變量（如t,u,v,w）等進(jìn)）等進(jìn)行說明，即要將這些變量說明成符號變量行說明，即要將這些變量說明成符號變量 Z變換例例.求數(shù)列求數(shù)列 fn=e-n的的Z變換及其逆變換。命令如下：變換及其逆變換。命令如下：syms n zfn=exp(-n);Fz=ztrans(fn,n,z) %求求fn的的Z變換變換f=iztrans(Fz,z,n) %求求Fz的逆的逆Z變換變換u例例 用用MATLAB求出離散序列求出離散序列 的的Z變換變換MATLAB程序如下：程序如下：syms k zf=0.5k; %定義離散信號定義離散信號Fz=ztrans(f) %對離散信號進(jìn)行對離散信號進(jìn)行Z變換變換u

9、運(yùn)行結(jié)果如下：運(yùn)行結(jié)果如下： Fz = 2*z/(2*z-1)Z變換u例例 已知一離散信號的已知一離散信號的Z變換式為變換式為 , 求出它所對應(yīng)的離散信號求出它所對應(yīng)的離散信號f(k). MATLAB程序如下：程序如下：syms k zFz=2*z/(2*z-1); %定義定義Z變換表達(dá)式變換表達(dá)式fk=iztrans(Fz,k) %求反求反Z變換變換u運(yùn)行結(jié)果如下：運(yùn)行結(jié)果如下：fk = (1/2)ku例例：求序列的：求序列的Z變換變換.syms n hn=sym(kroneckerDelta(n, 1) + kroneckerDelta(n, 2)+ kroneckerDelta(n, 3

10、)Hz=ztrans(hn)Hz=simplify(Hz)u運(yùn)行結(jié)果如下：運(yùn)行結(jié)果如下：Fz= (z2 + z + 1)/z3( )(1)(4)f kkt2( )21zF zz 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)的。 調(diào)用格式是: y=filter(b, a, x) x, y, a 和 b都是向量。與逆與逆Z變換變換 相關(guān)的相關(guān)的matlab 函數(shù)函數(shù)與逆與逆Z變換變換 相關(guān)的相

11、關(guān)的matlab 函數(shù)函數(shù)2.impz.m在 A（z）、B（z）已知情況下, 求系統(tǒng)的單位抽樣響應(yīng) h(n)。調(diào)用格式是： h = impz(b, a, N) 或 h,t=impz(b,a,N) N是所需的的長度。前者繪圖時(shí)n從1開始，而后者從0開始。 3. residuez.m將X(z) 的有理分式分解成簡單有理分式的和，因此可用來求逆變換。調(diào)用格式： r,p,k= residuez(b,a)假如知道了向量r, p和k，利用residuez.m還可反過來求出多項(xiàng)式A(z)、B(z)。格式是 b,a= residuez(r,p,k)。 頻率響應(yīng)函數(shù)：頻率響應(yīng)函數(shù)：freqz.m已知A(z)、B

12、(z), 求系統(tǒng)的頻率響應(yīng)?；镜恼{(diào)用格式是： H,w=freqz(b,a,N,whole,Fs)N是頻率軸的分點(diǎn)數(shù)，建議N為2的整次冪；w是返回頻率軸座標(biāo)向量，繪圖用；Fs是抽樣頻率，若Fs1，頻率軸給出歸一化頻率；whole指定計(jì)算的頻率范圍是從0FS,缺省時(shí)是從0FS/2.( )( )( )B zH zA z幅頻率響應(yīng)：幅頻率響應(yīng)：Hr=abs(H);相頻響應(yīng)：相頻響應(yīng)： Hphase=angle(H);解卷繞：解卷繞：Hphase=unwrap(Hphase);下面幾個(gè)文件用于轉(zhuǎn)移函數(shù)與極零點(diǎn)之 間的相互轉(zhuǎn)換及極零點(diǎn)的排序：（1) tf2zpk.m, 調(diào)用z,p,k=tf2zpk(b,

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

14、artImaginary Part1221222( )1 44(44)/(2) /H zzzzzzzz 用極零分析大致畫出幅頻響應(yīng)及相頻響應(yīng)：( )( )4 (1)4 (2)y nx nx nx n例1： 系統(tǒng)解：轉(zhuǎn)移函數(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;1221222( )1 44(44)/(2) /H zzzzzzzz 用極零分析大致畫出幅頻響應(yīng)及相頻響應(yīng)：( )( )4 (1

15、)4 (2)y nx nx nx n例1： 系統(tǒng)解：轉(zhuǎn)移函數(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

16、Res.|()| 1jH e11( )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-202405010015020025030002468-1-2-3-4-1-2-3-4.1836+.7344z +1.1016z +.7374z +.183116z( )1-3.0544z +3.8291z

17、-2.2925z +.55075z00H z 例：例： 給定系統(tǒng)給定系統(tǒng)求：極零圖 單位抽樣響應(yīng) 頻率響應(yī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解：-1.5-1-0.500.51-1-0.8-0.6-0.4-0.200.20.40.60.81

18、Real PartImaginary Part極零圖極零圖 zplane(b,a); 0510152025303540-0.1-0.0500.050.10.150.20.25單位抽樣響應(yīng)h,t=impz(b,a,40);stem(t,h,.);grid on;頻率響應(yīng)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

19、=angle(H); 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)下面幾個(gè)文件用于轉(zhuǎn)移函數(shù)、極零點(diǎn)與二階子系統(tǒng)二階子系統(tǒng)sos(Second-Order Section)級聯(lián)之間的相互轉(zhuǎn)換：（1) tf2sos.m, 調(diào)用sos,G=tf2sos(b,a) (2) sos2tf.m, 調(diào)用 b,a=sos2tf (sos,G) (3) sos2zp.m, 調(diào)用z,p,k=

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

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

22、、lp2hp.m、lp2bp.m, lp2bs.m將模擬低通原型轉(zhuǎn)換為實(shí)際的低通、高通、帶通及帶阻濾波器。b, a 是AF LP 的分子、分母的系數(shù)向量，B, A是轉(zhuǎn)換后的的分子、分母的系數(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是其帶寬。4bilinear.m ：雙線性變換，由模擬濾波器 得到數(shù)字濾波器。 Bz, Az=bilinear(B, A

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

24、) 用來設(shè)計(jì)數(shù)字帶阻濾波器，顯然，這時(shí)的Wn是12的向量；格式(4) 用來設(shè)計(jì)模擬濾波器。 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)例6.7.1(例6.5.1)uclear all; fp=100;fs=300;Fs=1000;rp=3;rs=20;uwp=2*pi*fp/Fs;ws=2*pi*fs/Fs;uFs=Fs/Fs; % let Fs=1 u% Firstly to finish frequency prewarping ;uwap=ta

25、n(wp/2);was=tan(ws/2); %un,wn=buttord(wap,was,rp,rs,s) % Note: s!uz,p,k=buttap(n); %ubp,ap=zp2tf(z,p,k) %ubs,as=lp2lp(bp,ap,wap) %u% Note: s=(2/Ts)(z-1)/(z+1);Ts=1,that is 2Fs=1,Fs=0.5;ubz,az=bilinear(bs,as,Fs/2) %uh,w=freqz(bz,az,256,Fs*1000);uplot(w,abs(h);grid on;設(shè)計(jì) IIR LP DF，3dB,20dBpssDF:100Hz,

26、300Hz,1000HzpsffF例6.7.1(例6.5.1)uclear all; uwp=.2*pi;ws=.6*pi;Fs=1000; rp=3;rs=20;%u% Firstly to finish frequency prewarping;uwap=2*Fs*tan(wp/2);was=2*Fs*tan(ws/2);un,wn=buttord(wap,was,rp,rs,s);% Note: s!uz,p,k=buttap(n);ubp,ap=zp2tf(z,p,k);ubs,as=lp2lp(bp,ap,wap)uw1=0:499*2*pi;uh1=freqs(bs,as,w1);

27、ubz,az=bilinear(bs,as,Fs) % Note: z=(2/ts)(z-1)/(z+1);uh2,w2=freqz(bz,az,500,Fs);uplot(w1/2/pi,abs(h1),w2,abs(h2),k);grid on;設(shè)計(jì) IIR LP DF，3dB,20dBpssDF:100Hz,300Hz,1000HzpsffF例6.7.1(例6.5.1)uclear all; uwp=.2*pi;ws=.6*pi;Fs=1000;urp=3;rs=20;un,wn=buttord(wp/pi,ws/pi,rp,rs);ubz,az=butter(n,wp/pi)ubz1,

28、az1=butter(n,wn)uh,w=freqz(bz,az,128,Fs);uh1,w1=freqz(bz1,az1,128,Fs);uplot(w,abs(h),w1,abs(h1),g.);grid on;設(shè)計(jì) IIR LP DF，3dB,20dBpssDF:100Hz,300Hz,1000HzpsffF非線性關(guān)系DF:100Hz,300Hz,1000Hz0.2 ,0.6 ,2AF:tan(/2)685.82109(Hz)2tan(/2)2452.762438(Hz)psspsppssssffFTT 設(shè)計(jì)的 AF 并不是按給定的技術(shù)指標(biāo)，但當(dāng)再由 變回 后，保證了 DF的技術(shù)要求。s

29、z又稱為頻率的預(yù)變形（Freq. Warping) 。例如 ： 2tan(/2)sT 抽樣頻率pf .2sf .2 給出給出數(shù)字高通數(shù)字高通的技術(shù)要求的技術(shù)要求)(zHdhpspsp, 得到得到模擬高通模擬高通的技術(shù)要求的技術(shù)要求)(sHahpspsp,)2tan(1step得到得到模擬低通模擬低通的技術(shù)要求的技術(shù)要求)(pGspsp,p12step設(shè)計(jì)出設(shè)計(jì)出)(pG3stepspp4step 得到得到模擬高通轉(zhuǎn)移模擬高通轉(zhuǎn)移函數(shù)函數(shù))(sHahp11zzs5step最后得到最后得到數(shù)字高通轉(zhuǎn)移數(shù)字高通轉(zhuǎn)移函數(shù)函數(shù))(zHdhp數(shù)字高通濾波器設(shè)計(jì)步驟7.6 數(shù)字高通, 帶通及帶阻濾波器的設(shè)計(jì)

30、psqpq/1對 帶通（BP）、帶阻（BS）數(shù)字濾波器的設(shè)計(jì)，只需改變圖中 Step2 和 Step4：22231BWBW 帶阻31213()sps 22231BWBW 帶通21331()sps 2N 121)(2pppG) 1() 1() 1(22222zzzpBW432142641. 0307. 1237. 2637. 11)21 (0201. 0)(zzzzzzzHdbps3dB,18dB,F2000Hzps要求：按上述轉(zhuǎn)換辦法，可以求出：例6.6.2： 設(shè)計(jì)一 IIR BP DF,要求：通帶頻率范圍 ： 300Hz 400Hz ;阻帶頻率范圍 ：200Hz、500Hz31 ffslsh

31、ff1buttord.m 確定 LP DF、或 LP AF的階次；（1） N, Wn = buttord(Wp, Ws, Rp, Rs)；（2）N, Wn = buttord(Wp, Ws, Rp, Rs,s)：與本章內(nèi)容有關(guān)的MATLAB文件(1)對應(yīng) 數(shù)字濾波器。其中 Wp, Ws分別是通帶和阻帶的截止頻率，其值在 01 之間，1對應(yīng)抽樣頻率的一半(歸一化頻率)。對低通和高通，Wp, Ws都是標(biāo)量，對帶通和帶阻，Wp, Ws是12的向量。Rp, Rs 分別是通帶和阻帶的衰減(dB)。N是求出的相應(yīng)低通濾波器的階次，Wn是求出的3dB頻率，它和Wp稍有不同。(2)對應(yīng)模擬濾波器，各變量含意和

32、（1）相同，但Wp, Ws及Wn的單位為弧度/秒，它們實(shí)際上是頻率。2buttap.m 設(shè)計(jì)模擬低通(Butt)原型濾波器。 z, p, k=buttap(N): N是欲設(shè)計(jì)的低通原型濾波器的階次，z, p, k是設(shè)計(jì)出的極點(diǎn)、零點(diǎn)及增益。3lp2lp.m、lp2hp.m、lp2bp.m, lp2bs.m將模擬低通原型轉(zhuǎn)換為實(shí)際的低通、高通、帶通及帶阻濾波器。b, a 是AF LP 的分子、分母的系數(shù)向量，B, A是轉(zhuǎn)換后的的分子、分母的系數(shù)向量；(1)中，Wo是低通或高通濾波器的截止頻率；B, A=lp2lp(b, a, Wo) B, A=lp2hp(b, a, Wo)(1)B, A=lp2

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

34、(1)(3) 用來設(shè)計(jì)數(shù)字濾波器，B，A分別是H(z)的分子、分母多項(xiàng)式的系數(shù)向量，Wn是通帶截止頻率，范圍在01之間。若Wn是標(biāo)量，(1)用來設(shè)計(jì)低通數(shù)字濾波器，若Wn是12的向量，則(1) 用來設(shè)計(jì)數(shù)字帶通濾波器； (2)用來設(shè)計(jì)數(shù)字高通濾波器； (3) 用來設(shè)計(jì)數(shù)字帶阻濾波器，顯然，這時(shí)的Wn是12的向量；格式(4) 用來設(shè)計(jì)模擬濾波器。 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)6cheb1ord.m 求Cheb-型濾波器的階；7cheb1a

35、p.m 設(shè)計(jì)原型低Cheb-I型模擬濾波器；8cheby1.m 直接設(shè)計(jì)數(shù)字Cheb-濾波器。以上三個(gè)文件的調(diào)用格式和對應(yīng)的Butterworth濾波器的文件類似。9cheb2ord.m; 10. ellipord.m; 11.cheb2ap.m; 12. ellipap.m; 13.besselap.m; 14. cheby2.m; 15. ellip.m; 16.besself.m17impinvar.m 用沖激響應(yīng)不變法實(shí)現(xiàn)頻率轉(zhuǎn)換； 對應(yīng) Cheby-II、橢圓 IIR 濾波器產(chǎn)生窗函數(shù)的文件有八個(gè)： bartlett（三角窗）; 2. blackman（布萊克曼窗） ; 3. box

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

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

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

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

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

41、,Wn, DC-0)前者保證第一個(gè)帶為通帶，后者保證第一個(gè)帶為阻帶。在上述所有格式中，若不指定窗函數(shù)的類型，fir1自動(dòng)選擇Hamming窗。指定窗函數(shù)格式： （4）b = fir1(N,Wn,wind); 例 b = fir1(N,Wn,boxcar(N+1); 指定矩形窗例7.1.1.設(shè)計(jì)低通 DF FIR, 令截止頻率 0. 25, 取 M10, 20，40，觀察其幅頻響應(yīng)的特點(diǎn).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,boxcar(4*N+1

42、); 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)()jdH e22cc1例7.1.2： 理想差分器及其設(shè)計(jì)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; b1(n)=0; e

43、lse 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 nnM11. remez

44、.m 設(shè)計(jì)Chebyshev最佳一致逼近FIR濾波器、Hilbert變換器和差分器。調(diào)用格式是： (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è)計(jì)的濾波器的系數(shù)，其長度為N1；F是頻率向量，A是對應(yīng)F的各頻段上的理想幅頻響應(yīng)，W是各頻段上的加權(quán)向量。F、A及W的指定方式和例7.4.1和7.4.2所討論過的一樣，唯一的差別是F的范圍為01，而非00.5, 1對應(yīng)抽樣頻率的一半。需要指出的

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

46、re(2);plot(w,h);grid; 調(diào)整通帶、阻帶的加權(quán)及濾波器的長度。調(diào)整N或W的結(jié)果例7.4.2： 設(shè)計(jì)多阻帶濾波器，抽樣頻率500Hz, 在50Hz、 100Hz 及150Hz處陷波。 通帶加權(quán)為8，阻帶為1-17dB通帶、阻帶加權(quán)都是1-25dB例7.4.2： 設(shè)計(jì)多阻帶濾波器，抽樣頻率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;weigh1=8 1 8 1 8 1 8;b1

47、=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.54 0.66250Hz12remezord.m 本文件用來確定在用Chebyshev最佳一致逼近設(shè)計(jì)FIR濾波器時(shí)所需要的濾波器階次。其調(diào)用格式是： N, Fo,

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

49、疊接相加法實(shí)現(xiàn)長序列卷積。格式是: y=fftfilt(h,x) 或 y=fftfilt(h, x，N)與本章有關(guān)的 MATLAB 文件記 的長度為 ， 的長度為 。 若采用第一個(gè)調(diào)用方式，程序自動(dòng)地確定對 分段的長度 及做FFT的長度 ， 顯然， 是最接近 的2的整次冪。分的段數(shù)為 。( )x nxN( )h nM( )x n()LMNL/xNLN采用第二個(gè)調(diào)用方式，使用者可自己指定做FFT的長度。建議使用第一個(gè)調(diào)用方式。clear;% 用疊接相加法，計(jì)算濾波器系數(shù)用疊接相加法，計(jì)算濾波器系數(shù)h和輸入信號和輸入信號x的卷積的卷積% 其中其中h為為10階階hanning窗，窗，x是帶有高斯白噪

50、的正弦信號是帶有高斯白噪的正弦信號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)plot(y);例3.9.1令x(n)為一正弦加

51、白噪聲信號，長度為500， h(n)是用fir1.m文件設(shè)計(jì)出的一個(gè)低通FIR濾波器，長度為11.試用fftfilt實(shí)現(xiàn)長序列的卷積實(shí)現(xiàn)長序列的卷積050100150200250300350400450500-2-1012050100150200250300350400450500-2-1012clear;% 生成濾波器系數(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:N)+s;% 將x分為長度為L的小段L=20;M=leng

52、th(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)050100150200250300350400450500-2-1012050100150200250300350400450500-2-1012 hilbert.m 文

53、件用來計(jì)算信號Hilbert變換。調(diào)用的格式是: y=hilbert(x),y的實(shí)部就是 ，虛部是的Hilbert變換 。所以，y 實(shí)際上是 x 的解析信號。( )x n( )x n czt.m 調(diào)用格式是： Xczt(x, M, W, A) 。x是待變換的時(shí)域信號，其長度設(shè)為N，M是變換的長度，W確定變換的步長，A確定變換的起點(diǎn)。若M=N, A=1, 則CZT變成DFT。00,jjAeWeA=exp(j*2*pi*f0/fs);W=exp(-j*2*pi*DELf/fs);例例。clear all; % 產(chǎn)生兩個(gè)正弦加白噪聲； N=256; f1=.1;f2=.2;fs=1; a1=5;a2

54、=3; w=2*pi/fs; x=a1*sin(w*f1*(0:N-1)+a2*sin(w*f2*(0:N-1)+randn(1,N); % 應(yīng)用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); 010203040506070-10010-0.5-0.4-0.3-0.2-0.100.10.20.30.40.502004006000102

55、03040506070-10010例例4.7.2 設(shè)設(shè)x(n)由三個(gè)實(shí)正弦組成，頻率分別由三個(gè)實(shí)正弦組成，頻率分別是是8HZ， 8.22HZ 和和9HZ， 抽樣頻率是抽樣頻率是40HZ ，時(shí)域時(shí)域取取128點(diǎn)。點(diǎn)。024681012141618200501000246810121416182005010077.588.599.51010.5050100 程序clear all；% 構(gòu)造三個(gè)不同頻率的正弦信號的疊加作為試驗(yàn)信號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/M);% A=1時(shí)的czt變換A=1;Y1=czt(x,M,W,A);subplot(311)plot(n1,abs(Y1(1:N/2);grid on;% DTFTY2=abs(fft(x);subplot(312)plot(n1,abs(Y2(1:N/2);grid on;% 詳細(xì)構(gòu)造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