Matlab與信號(hào)處理(信號(hào)處理)課件(PPT 80頁(yè))

1、MATLAB與信號(hào)處理數(shù)字信號(hào)處理.1第1頁(yè)，共80頁(yè)。4、基于MATLAB的數(shù)字信號(hào)處理4.1 信號(hào)產(chǎn)生4.2 離散系統(tǒng)的MATLAB實(shí)現(xiàn)4.3 變換的MATLAB實(shí)現(xiàn)4.4 數(shù)字濾波器設(shè)計(jì)4.5 功率譜估計(jì).2第2頁(yè)，共80頁(yè)。4.1 信號(hào)產(chǎn)生信號(hào)的表示x(n)=x(n)=,x(-2),x(-1),x(0),x(1),x(2) Matlab中，信號(hào)都是用矩陣（向量）描述例如：x=4 3 7 9 1典型離散信號(hào)表示單位抽樣序列 x=zeros(1,N); x(1)=1; (n) x(k)=1; (n-k)單位階躍序列: x=ones(1,N);.3第3頁(yè)，共80頁(yè)。波形的產(chǎn)生首先產(chǎn)生時(shí)間：如

2、t=0:0.01:2*pi;產(chǎn)生信號(hào)：如 y=sin(2*t+pi/4);square sawtooth sinc 噪聲的產(chǎn)生rand , randn.4第4頁(yè)，共80頁(yè)。 周期波形產(chǎn)生函數(shù)名功能Sawtooth(t,width)產(chǎn)生鋸齒波或三角波. Width確定最大值位置Square(t,duty)產(chǎn)生方波。Duty指定正半周期的比例Sinc(t)產(chǎn)生sinc或sin(pi*t)/(pi*t)Diric(t,n)產(chǎn)生Dirichlet函數(shù).5第5頁(yè)，共80頁(yè)。sawtooth(t,0.7)square(t)sinc(t)diric(t,5).6第6頁(yè)，共80頁(yè)。信號(hào)的常規(guī)操作信號(hào)加x=x1

3、+x2信號(hào)必須等長(zhǎng)對(duì)齊信號(hào)乘x=x1.*x2信號(hào)必須等長(zhǎng)對(duì)齊改變比例x=a*x1移位x=x1(n-k).7第7頁(yè)，共80頁(yè)。折疊y=x(-n)取樣和y=sum(x(n1:n2)取樣積y=prod(x(n1:n2)卷積y=conv(x,h)y,ny=conv_m(x,nx,h,nh)相關(guān)r= xcorr(x,y).8第8頁(yè)，共80頁(yè)。產(chǎn)生數(shù)據(jù)直接鍵盤輸入利用函數(shù)產(chǎn)生，如sin,cos,sawtooh等用Matlab從.Mat文件導(dǎo)入通過(guò)文件讀寫(xiě)產(chǎn)生數(shù)據(jù)，如fopen,fread等用一個(gè).M文件專門產(chǎn)生數(shù)據(jù).9第9頁(yè)，共80頁(yè)。4.2 離散系統(tǒng)的MATLAB實(shí)現(xiàn)離散系統(tǒng)的表示方法： LSI系統(tǒng)的時(shí)

4、域表示 常系數(shù)線性差分方程 單位抽樣響應(yīng)h(n) LSI系統(tǒng)的頻域表示 頻率響應(yīng) （DFT） 轉(zhuǎn)移函數(shù) （z變換） 零極點(diǎn)增益 二次分式 離散系統(tǒng)的內(nèi)部描述 .10第10頁(yè)，共80頁(yè)。離散系統(tǒng)的輸入輸出的差分方程：?jiǎn)挝怀闃禹憫?yīng) ?filter函數(shù)y=filter(b,a,x)impz函數(shù)y=impz(b,a).11第11頁(yè)，共80頁(yè)。.12第12頁(yè)，共80頁(yè)。頻率響應(yīng) ？ h,f=freqz(b,a,n,fs).13第13頁(yè)，共80頁(yè)。.14第14頁(yè)，共80頁(yè)。 零極點(diǎn)增益？.15第15頁(yè)，共80頁(yè)。離散系統(tǒng)的表示方法轉(zhuǎn)換函數(shù) tf2zp sos2tf tf2ss sos2zp zp2tf s

5、os2ss zp2sos ss2tf zp2ss ss2zp ss2sos.16第16頁(yè)，共80頁(yè)。 4.3 變換的MATLAB實(shí)現(xiàn)1、DFT.17第17頁(yè)，共80頁(yè)。計(jì)算 DFT并畫(huà)圖 .18第18頁(yè)，共80頁(yè)。.19第19頁(yè)，共80頁(yè)。.20第20頁(yè)，共80頁(yè)。2、FFT 可利用內(nèi)部函數(shù)fft進(jìn)行計(jì)算，速度快。 y=fft(x,n) x的n點(diǎn)fft x=ifft(y,n).21第21頁(yè)，共80頁(yè)。3、z變換R,P,K=residuez(b,a);.22第22頁(yè)，共80頁(yè)。4、DCTy=dct(x,n)x=idct(y,n).23第23頁(yè)，共80頁(yè)。rp=84.3566.24第24頁(yè)，共80

6、頁(yè)。5、Hilbert變換 解析信號(hào) DFT方法求解 函數(shù)y=Hilbert(x).25第25頁(yè)，共80頁(yè)。（1） 數(shù)字濾波器的基本概念（2）基于MATLAB的IIR濾波器設(shè)計(jì)（3） 基于MATLAB的FIR濾波器設(shè)計(jì)4.4 基于MATLAB的數(shù)字濾波器設(shè)計(jì).26第26頁(yè)，共80頁(yè)。（1） 數(shù)字濾波器的基本概念數(shù)字濾波器 對(duì)數(shù)字信號(hào)進(jìn)行濾波的線性時(shí)不變系統(tǒng)。數(shù)字濾波器的分類 IIR數(shù)字濾波器：設(shè)計(jì)經(jīng)典選頻濾波器，性能好， 階數(shù)低。 FIR數(shù)字濾波器：設(shè)計(jì)線性相位經(jīng)典濾波器，適用范圍廣、 階數(shù)高。 .27第27頁(yè)，共80頁(yè)。數(shù)字濾波器的設(shè)計(jì)步驟1） 按照實(shí)際任務(wù)的要求，確定濾波器的技術(shù)指標(biāo)。2）

7、用一個(gè)因果、穩(wěn)定的離散線性時(shí)不變系統(tǒng)的系統(tǒng)函數(shù) 去逼近這一網(wǎng)絡(luò)性能。根據(jù)不同的要求可選擇IIR或者 FIR系統(tǒng)函數(shù)去逼近。3）利用有限精度算法實(shí)現(xiàn)系統(tǒng)函數(shù)。.28第28頁(yè)，共80頁(yè)。IIR數(shù)字濾波器一般有兩種設(shè)計(jì)方法： 1）模擬原型法 先設(shè)計(jì)一個(gè)合適的模擬濾波器，然后再變 換為滿足要求的數(shù)字濾波器 。 2）直接法 直接在離散時(shí)域和離散頻域進(jìn)行數(shù)字濾波 器的設(shè)計(jì) .29第29頁(yè)，共80頁(yè)。FIR一般有三種設(shè)計(jì)方法： 窗函數(shù)法 頻率取樣法 切比雪夫逼近法 .30第30頁(yè)，共80頁(yè)。（2） 基于MATLAB的IIR濾波器設(shè)計(jì)Butterworth數(shù)字濾波器設(shè)計(jì)橢圓數(shù)字濾波器設(shè)計(jì)Chebyshev-I

8、數(shù)字濾波器設(shè)計(jì)Chebyshev-II數(shù)字濾波器設(shè)計(jì).31第31頁(yè)，共80頁(yè)。Butterworth數(shù)字濾波器設(shè)計(jì) MATLAB語(yǔ)句： b,a=butter(n,wn) (1) b,a=butter(n,wn,ftype) (2) z,p,k=butter() (3) a,b,c,d=butter() (4) 語(yǔ)句(1) (4) 都可以設(shè)計(jì)一個(gè)階數(shù)為n，截止頻率為wn的數(shù)字低通濾波器。對(duì)于語(yǔ)句(2) ，若參數(shù)ftype取hi-gh或stop，則可以來(lái)設(shè)計(jì)相應(yīng)的高通或帶阻濾波器。 返回值a和b為系統(tǒng)函數(shù)的分子和分母的系數(shù)。.32第32頁(yè)，共80頁(yè)。系統(tǒng)函數(shù)為： 返回值z(mì),p,k分別是濾波器的零點(diǎn)

9、、極點(diǎn)和增益。返回值A(chǔ)、B、C、D構(gòu)造濾波器的狀態(tài)方程。狀態(tài)方程為： 其中u是輸入信號(hào)，x是狀態(tài)變量，y是輸出信號(hào)。.33第33頁(yè)，共80頁(yè)。例 ： 設(shè)計(jì)一個(gè)階數(shù)為9階、截止頻率為300Hz的高通 Butterworth濾波器，采樣頻率為1000。并畫(huà)出它的頻率特性。程序： n=9; wn=300/500; b,a=buffer(n,wn,high); freqz(b,a,128,1000);.34第34頁(yè)，共80頁(yè)。Buttordworth高通濾波器的頻率響應(yīng).35第35頁(yè)，共80頁(yè)。 n=2.36第36頁(yè)，共80頁(yè)。例 設(shè)計(jì)一個(gè)10階的帶通Butterworth濾波器，它的通帶 范圍是10

10、0200Hz,采樣頻率為1000。并畫(huà)出它的頻率特性和沖擊響應(yīng)。 程序：n=10; wn=100 200/500; b,a=butter(n,wn); freqz(b,a,128,1000); y,t=impz(b,a,101); stem(t,y);.37第37頁(yè)，共80頁(yè)。 Buttordworth帶通濾波器的頻率響應(yīng).38第38頁(yè)，共80頁(yè)。Buttordworth濾波器的時(shí)域沖擊響應(yīng).39第39頁(yè)，共80頁(yè)。橢圓數(shù)字濾波器設(shè)計(jì) MATLAB語(yǔ)法： b,a=ellip(n,Rp,Rs,wn) (5) b,a=ellip(n, Rp,Rs, wn,ftype) (6) z,p,k=elli

11、p() (7) a,b,c,d=ellip() (8) 語(yǔ)句(5) (8)都可以設(shè)計(jì)一個(gè)階數(shù)為n，截止頻率為wn,通帶波紋最大衰減為Rp、阻帶波紋最小衰減為Rs的數(shù)字低通濾波器。對(duì)于語(yǔ)句(6)，若參數(shù)ftype取high或stop，則可以來(lái)設(shè)計(jì)相應(yīng)的高通或帶阻濾波器。 返回值a和b為系統(tǒng)函數(shù)的分子和分母的系數(shù)。返回 .40第40頁(yè)，共80頁(yè)。值z(mì),p,k分別是濾波器的零點(diǎn)、極點(diǎn)和增益。返回值A(chǔ)、B、C、D可以用來(lái)構(gòu)造濾波器的狀態(tài)方程。例 采樣頻率為1000Hz，請(qǐng)?jiān)O(shè)計(jì)一個(gè)階數(shù)為6階、截止頻 率為300Hz的低通ellip數(shù)字濾波器。其中濾波器在通 帶的波紋為3dB，阻帶的波紋為50dB。 程序

12、： b,a=ellip(6,3,50,300/500); freqz(b,a,512,1000); title(n=6 低通橢圓濾波器);.41第41頁(yè)，共80頁(yè)。橢圓低通濾波器的頻率響應(yīng).42第42頁(yè)，共80頁(yè)。例 設(shè)計(jì)一個(gè)10階的帶通濾波器，通帶范圍為100200Hz， 信號(hào)的采樣頻率為1000Hz，濾波器在通帶的波紋為0.5 dB，阻帶的波紋為20dB。設(shè)計(jì)出這個(gè)濾波器并畫(huà)出它 的沖擊響應(yīng)。 程序： b,a=ellip(10,0.5,20,100 200/500); y,t=impz(b,a,101); freqz(b,a,128,1000); stem(t,y); title(n=10

13、 帶通橢圓濾波器);.43第43頁(yè)，共80頁(yè)。橢圓帶通濾波器的頻率響應(yīng).44第44頁(yè)，共80頁(yè)。橢圓濾波器的沖擊響應(yīng).45第45頁(yè)，共80頁(yè)。Chebyshev-I數(shù)字濾波器設(shè)計(jì) MATLAB語(yǔ)法： b,a=cheby1(n,Rp,wn) (9) b,a= cheby1(n,Rp,wn,ftype) (10) z,p,k= cheby1() (11) a,b,c,d= cheby1() (12) 語(yǔ)句(9) (12)都可以設(shè)計(jì)一個(gè)階數(shù)為n，截止頻率為wn,通帶波紋最大衰減為Rp的數(shù)字低通濾波器。對(duì)于語(yǔ)句(10)，若參數(shù)ftype取high或stop，則可以來(lái)設(shè)計(jì)相應(yīng)的高通或帶阻濾波器。 返回值

14、a和b為系統(tǒng)函數(shù)的分子和分母的系數(shù)。返回.46第46頁(yè)，共80頁(yè)。值z(mì),p,k分別是濾波器的零點(diǎn)、極點(diǎn)和增益。返回值A(chǔ)、B、C、D可以用來(lái)構(gòu)造濾波器的狀態(tài)方程。例 對(duì)采樣頻率為1000Hz的采樣信號(hào)，設(shè)計(jì)一個(gè)階數(shù)為9 階，截止頻率為300Hz的低通cheby1數(shù)字濾波器，其 中濾波器在通帶的波紋為0.5dB。 程序： b,a=cheby1(9,0.5,300/500); freqz(b,a,512,1000); .47第47頁(yè)，共80頁(yè)。 Chebyshev-I數(shù)字低通濾波器的頻率響應(yīng).48第48頁(yè)，共80頁(yè)。Chebyshev-II數(shù)字濾波器設(shè)計(jì) MATLAB語(yǔ)法： b,a=cheby2(n

15、,Rs,wn) (13) b,a= cheby2(n,Rs,wn,ftype) (14) z,p,k= cheby2() (15) a,b,c,d= cheby2() (16) 語(yǔ)句(13) (16)都可以設(shè)計(jì)一個(gè)階數(shù)為n，截止頻率為wn, 阻帶紋波最小衰減為Rs的數(shù)字低通濾波器。對(duì)于語(yǔ)句(13)式，若參數(shù)ftype取high或stop，則可以來(lái)設(shè)計(jì)相應(yīng)的高通或帶阻濾波器。 返回值a和b為系統(tǒng)函數(shù)的分子和分母的系數(shù)。返回.49第49頁(yè)，共80頁(yè)。值z(mì),p,k分別是濾波器的零點(diǎn)、極點(diǎn)和增益。返回值A(chǔ)、B、C、D可以用來(lái)構(gòu)造濾波器的狀態(tài)方程。例 對(duì)于采樣頻率為1000Hz的信號(hào)，設(shè)計(jì)一個(gè)階數(shù)為9階

16、，截止頻率為300Hz的低通Chebyshev-II數(shù)字濾波器，其中濾波器在阻帶的波紋為20dB。 程序： b,a=cheby2(9,20,300/500); freqz(b,a,512,1000); .50第50頁(yè)，共80頁(yè)。 Chebyshev-II數(shù)字低通濾波器的頻率響應(yīng) .51第51頁(yè)，共80頁(yè)。（3 ）基于MATLAB的FIR濾波器設(shè)計(jì)窗函數(shù)法設(shè)計(jì)FIR數(shù)字濾波器頻率取樣法切比雪夫逼近法設(shè)計(jì)FIR濾波器.52第52頁(yè)，共80頁(yè)。 窗函數(shù)法設(shè)計(jì)FIR數(shù)字濾波器 時(shí)域中進(jìn)行的 常見(jiàn)的窗函數(shù) w=boxcar(n) w=triang(n) w=hanning(n) w=hamming(n)

17、 w=blackman(n) w=chebwin(n,R) w=bartlett(n) w=kaiser(n,) .53第53頁(yè)，共80頁(yè)。例：用8種窗函數(shù)設(shè)計(jì)數(shù)字低通濾波器。采樣頻率 1000Hz，截止頻率200Hz，濾波器階數(shù)81。 程序：.54第54頁(yè)，共80頁(yè)。.55第55頁(yè)，共80頁(yè)。利用fir1函數(shù)設(shè)計(jì)濾波器（窗函數(shù)法） MATLAB語(yǔ)法： b=fir1(n, wn) b=fir1(n, wn, ftype) b=fir1(n, wn, window) b=fir1(n, wn, ftype, window) .56第56頁(yè)，共80頁(yè)。例 : 設(shè)計(jì)一個(gè)階數(shù)為48，通帶范圍為 0.3

18、50.65的帶通FIR線性相位濾波器，并分析它的頻率特性。 程序：b=fir1(48,0.35,0.65); freqz(b); .57第57頁(yè)，共80頁(yè)。 利用FIR1設(shè)計(jì)的帶通濾波器的頻率響應(yīng).58第58頁(yè)，共80頁(yè)。利用kaiserord函數(shù)求凱賽窗函數(shù)的參數(shù) MATLAB語(yǔ)法： n,wn,beta,ftype=kaiserord(f,a,dev) n,wn,beta,ftype=kaiserord(f,a,dev,Fs）例 利用凱賽窗函數(shù)設(shè)計(jì)一個(gè)低通FIR數(shù)字濾波器，通帶 范圍是01000Hz，阻帶范圍是15004000Hz，通帶 的波紋最大為0.05，阻帶的波紋最大為0.01。采樣頻

19、 率為8000Hz。.59第59頁(yè)，共80頁(yè)。 程序: fsamp=8000; fcuts=1000 1500; mags=1 0; devs=0.05 0.01; n,wn,beta,ftype=kaiserord(fcuts,mags, devs,fsamp); b=fir1(n,wn,ftype,kaiser(n+1,beta); freqz(b); .60第60頁(yè)，共80頁(yè)。 利用凱賽窗設(shè)計(jì)的低通濾波器的頻率響應(yīng).61第61頁(yè)，共80頁(yè)。利用fir2函數(shù)設(shè)計(jì)任意響應(yīng)FIR數(shù)字濾波器 MATLAB語(yǔ)法： b=fir2(n, f, m) b=fir2(n, f, m, window) b=

20、fir2(n, f, m, npt) b=fir2(n, f, m, npt, window) b=fir2(n, f, m, npt, lap) .62第62頁(yè)，共80頁(yè)。例 設(shè)計(jì)一個(gè)60階的濾波器，要求設(shè)計(jì)的濾波器在0 的幅度響應(yīng)為1，在 到 的幅度響應(yīng)為 ，在 到 的幅度響應(yīng)為 ，在 到 的 幅度響應(yīng)為 ，在 到 的幅度響應(yīng)為 。并畫(huà) 出該濾波器和理想濾波器的幅頻響應(yīng)曲線，進(jìn)行比較。 程序：f=0,0.125,0.125,0.25,0.25,0.5,0.5,0.75,0.75,1; m=1,1,0.5,0.5,0.25,0.25,1/6,1/6,0.125,0.125; b=fir2(6

21、0, f, m); h,w= freqz(b); plot(f,m,w/pi,abs(h);.63第63頁(yè)，共80頁(yè)。 用fir2設(shè)計(jì)的任意響應(yīng)的濾波器與理想濾波器比較.64第64頁(yè)，共80頁(yè)。 切比雪夫逼近法設(shè)計(jì)FIR濾波器 采用最大誤差最小準(zhǔn)則方法設(shè)計(jì)濾波器 優(yōu)化問(wèn)題 雷米茲（REMEZ)算法 較成功的設(shè)計(jì)線性相位濾波器的方法.65第65頁(yè)，共80頁(yè)。利用remez函數(shù)設(shè)計(jì)FIR數(shù)字濾波器 MATLAB語(yǔ)法： b=remez(n, f, m) b=remez(n, f, m, w) b=remez(n, f, m, ftype) b=remez(n, f, m, w, ftype) .66

22、第66頁(yè)，共80頁(yè)。其它的MATLAB函數(shù) ： firls remezord fircs fircls1 cremez firrcos .67第67頁(yè)，共80頁(yè)。 IIR與FIR數(shù)字濾波器的比較在相同的技術(shù)指標(biāo)下，IIR可以用比FIR較少的階數(shù)來(lái)滿足要求；FIR可得到嚴(yán)格的線性相位；IIR的設(shè)計(jì)可利用模擬濾波器的結(jié)果IIR主要設(shè)計(jì)規(guī)格化的濾波器，而FIR較靈活，可是任意響應(yīng)的濾波器，適應(yīng)性較廣。 .68第68頁(yè)，共80頁(yè)。一維信號(hào)處理：信號(hào)濾波信號(hào)：s1 = sin(2*pi*t*5); s2=sin(2*pi*t*15); s3=sin(2*pi*t*30);S=s1+s2+s3濾波器：b,

23、a = ellip(4,0.1,40,10 20*2/Fs);H,w = freqz(b,a,512);濾波后的結(jié)果sf = filter(b,a,s);.69第69頁(yè)，共80頁(yè)。5 功率譜估計(jì) 概念 ： 平穩(wěn)隨機(jī)信號(hào)的功率譜密度（PSD）是自相關(guān)序列的離散時(shí)間傅立葉變換.70第70頁(yè)，共80頁(yè)。 功率譜估計(jì)的方法.71第71頁(yè)，共80頁(yè)。經(jīng)典功率譜估計(jì)方法直接法（周期圖法） 直接由傅立葉變換而得來(lái)的。.72第72頁(yè)，共80頁(yè)。例子：直接法計(jì)算功率譜利用FFT直接法：Fs=1000;NFFT=1024;n=0:1/Fs:1;x=sin(2*pi*40*n)+4*sin(2*pi*100*n)+randn(size(n); %噪聲序列X=fft(x,NFFT);Pxx=abs(X).2/length(n); %求解PSDt=0:round(N