.16/161數(shù)字信號處理各章節(jié)重點知識的matlab實例第1,2章離散時間信號與系統(tǒng)例1-1用MATLAB計算序列{-2

0

1

–1

3}和序列{1

2

0

-1}的離散卷積。解MATLAB程序如下：

a=[-201-13];

b=[120-1];

c=conv<a,b>;

M=length<c>-1;

n=0:1:M;stem<n,c>;xlabel<'n'>;ylabel<'幅度'>;圖1.1給出了卷積結果的圖形,求得的結果存放在數(shù)組c中為：{-2

-4

1

3

1

5

1

-3}。例1-2用MATLAB計算差分方程當輸入序列為時的輸出結果。解MATLAB程序如下：

N=41;

a=[0.8-0.440.360.22];

b=[10.7-0.45-0.6];

x=[1zeros<1,N-1>];

k=0:1:N-1;

y=filter<a,b,x>;stem<k,y>

xlabel<'n'>;ylabel<'幅度'>圖1.2給出了該差分方程的前41個樣點的輸出,即該系統(tǒng)的單位脈沖響應。例1-3用MATLAB計算例1-2差分方程所對應的系統(tǒng)函數(shù)的DTFT。解例1-2差分方程所對應的系統(tǒng)函數(shù)為：其DTFT為

用MATLAB計算的程序如下：

k=256;

num=[0.8-0.440.360.02];

den=[10.7-0.45-0.6];

w=0:pi/k:pi;

h=freqz<num,den,w>;subplot<2,2,1>;

plot<w/pi,real<h>>;grid

title<'實部'>

xlabel<'\omega/\pi'>;ylabel<'幅度'>

subplot<2,2,2>;

plot<w/pi,imag<h>>;grid

title<'虛部'>

xlabel<'\omega/\pi'>;ylabel<'Amplitude'>

subplot<2,2,3>;

plot<w/pi,abs<h>>;grid

title<'幅度譜'>

xlabel<'\omega/\pi'>;ylabel<'幅值'>

subplot<2,2,4>;

plot<w/pi,angle<h>>;grid

title<'相位譜'>

xlabel<'\omega/\pi'>;ylabel<'弧度'>第3、4章離散傅里葉變換及其快速算法例2-1對連續(xù)的單一頻率周期信號,按采樣頻率采樣,截取長度N分別選N=20和N=16,觀察其DFT結果的幅度譜。解此時離散序列,即k=8。用MATLAB計算并作圖,函數(shù)fft用于計算離散傅里葉變換DFT,程序如下：

k=8;

n1=[0:1:19];

xa1=sin<2*pi*n1/k>;

subplot<2,2,1>

plot<n1,xa1>

xlabel<'t/T'>;ylabel<'x<n>'>;

xk1=fft<xa1>;xk1=abs<xk1>;

subplot<2,2,2>

stem<n1,xk1>

xlabel<'k'>;ylabel<'X<k>'>;

n2=[0:1:15];

xa2=sin<2*pi*n2/k>;

subplot<2,2,3>

plot<n2,xa2>

xlabel<'t/T'>;ylabel<'x<n>'>;

xk2=fft<xa2>;xk2=abs<xk2>;

subplot<2,2,4>

stem<n2,xk2>

xlabel<'k'>;ylabel<'X<k>'>;計算結果示于圖2.1,<a>和<b>分別是N=20時的截取信號和DFT結果,由于截取了兩個半周期,頻譜出現(xiàn)泄漏；<c>和<d>分別是N=16時的截取信號和DFT結果,由于截取了兩個整周期,得到單一譜線的頻譜。上述頻譜的誤差主要是由于時域中對信號的非整周期截斷產(chǎn)生的頻譜泄漏。例2-2用FFT計算兩個序列的互相關函數(shù)。解用MATLAB計算程序如下：

x=[13-112331];

y=[21-1120-13];

k=length<x>;

xk=fft<x,2*k>;

yk=fft<y,2*k>;

rm=real<ifft<conj<xk>.*yk>>;

rm=[rm<k+2:2*k>rm<1:k>];

m=<-k+1>:<k-1>;stem<m,rm>

xlabel<'m'>;ylabel<'幅度'>;其計算結果如圖2.2所示。例2-3計算兩個序列的的互相關函數(shù),其中x<n>={23521–100123530–1–2012}y<n>=x<n-4>+e<n>,e<n>為一隨機噪聲,在MATLAB中可以用隨機函數(shù)rand產(chǎn)生解用MATLAB計算程序如下：

x=[23521-100123530-1-2012];

y=[000023521-100123530-1-2012];

k=length<y>;

e=rand<1,k>-0.5;

y=y+e;

xk=fft<x,2*k>;

yk=fft<y,2*k>;

rm=real<ifft<conj<xk>.*yk>>;

rm=[rm<k+2:2*k>rm<1:k>];

m=<-k+1>:<k-1>;stem<m,rm>

xlabel<'m'>;ylabel<'幅度'>;計算結果如圖2.3<a>,我們看到最大值出現(xiàn)在m=4處,正好是y<n>對于x<n>的延遲。2.3<b>是x<n>的自相關函數(shù),他和y<n>的區(qū)別除時間位置外,形狀也略不同,這是由于y<n>受到噪聲的干擾。第5章無限長單位脈沖響應<IIR>濾波器的設計方法例3-1設采樣周期T=250μs〔采樣頻率fs=4kHz,用脈沖響應不變法和雙線性變換法設計一個三階巴特沃茲濾波器,其3dB邊界頻率為fc=1kHz。 [B,A]=butter<3,2*pi*1000,'s'>; [num1,den1]=impinvar<B,A,4000>; [h1,w]=freqz<num1,den1>; [B,A]=butter<3,2/0.00025,'s'>; [num2,den2]=bilinear<B,A,4000>; [h2,w]=freqz<num2,den2>; f=w/pi*2000; plot<f,abs<h1>,'-.',f,abs<h2>,'-'>; grid; xlabel<'頻率/Hz'> ylabel<'幅值/dB'>程序中第一個butter的邊界頻率2π×1000,為脈沖響應不變法原型低通濾波器的邊界頻率；第二個butter的邊界頻率2/T=2/0.00025,為雙線性變換法原型低通濾波器的邊界頻率.圖3.1給出了這兩種設計方法所得到的頻響,虛線為脈沖響應不變法的結果；實線為雙線性變換法的結果。脈沖響應不變法由于混疊效應,使得過渡帶和阻帶的衰減特性變差,并且不存在傳輸零點。同時,也看到雙線性變換法,在z=-1即ω=π或f=2000Hz處有一個三階傳輸零點,這個三階零點正是模擬濾波器在Ω=∞處的三階傳輸零點通過映射形成的。例3-2設計一數(shù)字高通濾波器,它的通帶為400～500Hz,通帶內(nèi)容許有0.5dB的波動,阻帶內(nèi)衰減在小于317Hz的頻帶內(nèi)至少為19dB,采樣頻率為1,000Hz。 wc=2*1000*tan<2*pi*400/<2*1000>>; wt=2*1000*tan<2*pi*317/<2*1000>>; [N,wn]=cheb1ord<wc,wt,0.5,19,'s'>; [B,A]=cheby1<N,0.5,wn,'high','s'>; [num,den]=bilinear<B,A,1000>; [h,w]=freqz<num,den>; f=w/pi*500; plot<f,20*log10<abs<h>>>; axis<[0,500,-80,10]>; grid; xlabel<''> ylabel<'幅度/dB'>圖3.2給出了MATLAB計算的結果,可以看到模擬濾波器在Ω=∞處的三階零點通過高通變換后出現(xiàn)在ω=0〔z=1處,這正是高通濾波器所希望得到的。例3-3設計一巴特沃茲帶通濾波器,其３dB邊界頻率分別為f2=110kHz和f1=90kHz,在阻帶f3=120kHz處的最小衰減大于１０dB,采樣頻率fs=400kHz。w1=2*400*tan<2*pi*90/<2*400>>; w2=2*400*tan<2*pi*110/<2*400>>; wr=2*400*tan<2*pi*120/<2*400>>; [N,wn]=buttord<[w1w2],[0wr],3,10,'s'>; [B,A]=butter<N,wn,'s'>; [num,den]=bilinear<B,A,400>; [h,w]=freqz<num,den>; f=w/pi*200; plot<f,20*log10<abs<h>>>; axis<[40,160,-30,10]>; grid; xlabel<'頻率/kHz'> ylabel<'幅度/dB'>圖3.3給出了MATLAB計算的結果,可以看出數(shù)字濾波器將無窮遠點的二階零點映射為z=±1的二階零點,數(shù)字帶通濾波器的極點數(shù)是模擬低通濾波器的極點數(shù)的兩倍。例3-4一數(shù)字濾波器采樣頻率fs=1kHz,要求濾除100Hz的干擾,其３dB的邊界頻率為95Hz和105Hz,原型歸一化低通濾波器為 w1=95/500; w2=105/500; [B,A]=butter<1,[w1,w2],'stop'>; [h,w]=freqz<B,A>; f=w/pi*500; plot<f,20*log10<abs<h>>>; axis<[50,150,-30,10]>; grid; xlabel<'頻率/Hz'> ylabel<'幅度/dB'>圖3.4為MATLAB的計算結果第5章有限長單位脈沖響應<FIR>濾波器的設計方法例4.1用凱塞窗設計一FIR低通濾波器,低通邊界頻率,阻帶邊界頻率,阻帶衰減不小于50dB。解首先由過渡帶寬和阻帶衰減來決定凱塞窗的N和圖4.1給出了以上設計的頻率特性,<a>為N=30直接截取的頻率特性<b>為凱塞窗設計的頻率特性。凱塞窗設計對應的MATLAB程序為： wn=kaiser<30,4.55>; nn=[0:1:29]; alfa=<30-1>/2; hd=sin<0.4*pi*<nn-alfa>>./<pi*<nn-alfa>>; h=hd.*wn'; [h1,w1]=freqz<h,1>; plot<w1/pi,20*log10<abs<h1>>>; axis<[0,1,-80,10]>; grid; xlabel<'歸一化頻率/'> ylabel<'幅度/dB'>例4-2利用雷米茲交替算法,設計一個線性相位低通FIR數(shù)字濾波器,其指標為：通帶邊界頻率fc=800Hz,阻帶邊界fr=1000Hz,通帶波動阻帶最小衰減At=40dB,采樣頻率fs=4000Hz。解在MATLAB中可以用remezord和remez兩個函數(shù)設計,其結果如圖4.2,MATLAB程序如下： fedge=[8001000]; mval=[10]; dev=[0.05590.01];

fs=4000; [N,fpts,mag,wt]=remezord<fedge,mval,dev,fs>; b=remez<N,fpts,mag,wt>; [h,w]=freqz<b,1,256>; plot<w*2000/pi,20*log10<abs<h>>>; grid;

xlabel<'頻率/Hz'> ylabel<'幅度/dB'>函數(shù)remezord中的數(shù)組fedge為通帶和阻帶邊界頻率,數(shù)組mval是兩個邊界處的幅值,而數(shù)組dev是通帶和阻帶的波動,fs是采樣頻率單位為Hz。6章數(shù)字信號處理系統(tǒng)的實現(xiàn)例5-1求下列直接型系統(tǒng)函數(shù)的零、極點,并將它轉(zhuǎn)換成二階節(jié)形式解用MATLAB計算程序如下： num=[1-0.1-0.3-0.3-0.2]; den=[10.10.20.20.5]; [z,p,k]=tf2zp<num,den>; m=abs<p>; disp<'零點'>;disp<z>; disp<'極點'>;disp<p>; disp<'增益系數(shù)'>;disp<k>; sos=zp2sos<z,p,k>; disp<'二階節(jié)'>;disp<real<sos>>; zplane<num,den>輸入到"num"和"den"的分別為分子和分母多項式的系數(shù)。計算求得零、極點增益系數(shù)和二階節(jié)的系數(shù)：零點0.9615-0.5730-0.1443+0.5850i-0.1443-0.5850i極點0.5276+0.6997i0.5276-0.6997i-0.5776+0.5635i-0.5776-0.5635i增益系數(shù)1二階節(jié)1.0000-0.3885-0.55091.00001.15520.65111.00000.28850.36301.0000-1.05520.7679系統(tǒng)函數(shù)的二階節(jié)形式為：極點圖見圖5.1。例5-2分析五階橢圓低通濾波器的量化效應,其截止頻率為0.4,通帶紋波為0.4dB,最小的阻帶衰減為50dB。對濾波器進行截尾處理時,使用函數(shù)a2dT.m.。解用以下MATLAB程序分析量化效應 clf; [b,a]=ellip<5,0.4,50,0.4>; [h,w]=freqz<b,a,512>; g=20*log10<abs<h>>; bq=a2dT<b,5>; aq=a2dT<a,5>; [hq,w]=freqz<bq,aq,512>;

gq=20*log10<abs<hq>>; plot<w/pi,g,'b',w/pi,gq,'r:'>; grid; axis<[01-805]>; xlabel<'\omega/\pi'>; ylabel<'Gain,dB'>; legend<'量化前','量化后'>; figure [z1,p1,k1]=tf2zp<b,a>; [z2,p2,k2]=tf2zp<bq,aq>; zplaneplot<[z1,z2],[p1,p2],{'o','x','*','+'}>; legend<'量化前的零點','量化后的零點','量化前的極點','量化后的極點'>;圖5.1〔a表示系數(shù)是無限精度的理想濾波器的頻率響應〔以實線表示以及當濾波器系數(shù)截尾到5位時的頻率響應〔以短線表示。由圖可知,系數(shù)量化對頻帶的邊緣影響較大,經(jīng)系數(shù)量化后,增加了通帶的波紋幅度,減小了過渡帶寬,并且減小了最小的阻帶衰減。圖5.1〔b給出了系數(shù)量化以前和系數(shù)量化以后的橢圓低通濾波器的零極點位置。由圖可知,系數(shù)的量化會使零極點的位置與它們的理想的標稱位置相比發(fā)生顯著的改變。在這個例子中,靠近虛軸的零點的位置變動最大,并且移向靠它最近的極點的位置。只要對程序稍作改變就可以分析舍入量化的影響。為了研究二進制數(shù)量化效應對數(shù)字濾波器的影響,首先需要將十進制表示的濾波器系數(shù)轉(zhuǎn)換成二進制數(shù)并進行量化,二進制數(shù)的量化既可以通過截尾法也可以通過舍入法實現(xiàn)。我們提供了如下的兩個MATLAB程序：a2dT.m和a2dR.m,這兩段程序分別將向量d中的每一個數(shù)按二進制數(shù)進行截尾或舍入量化,量化的精度是小數(shù)點以后保留b位,量化后返回的向量為beq。functionbeq=a2dT<d,b>%beq=a2dT<d,b>將十進制數(shù)利用截尾法得到b位的二進制數(shù),%然后將該二進制數(shù)再轉(zhuǎn)換為十進制數(shù)m=1;d1=abs<d>;whilefix<d1>>0d1=abs<d>/<2^m>;m=m+1;endbeq=fix<d1*2^b>;beq=sign<d>.*beq.*2^<m-b-1>;functionbeq=a2dR<d,b>%beq=a2dR<d,b>將十進制數(shù)利用舍入法得到b位的二進制數(shù)%然后將該二進制數(shù)再轉(zhuǎn)換為十進制數(shù)m=1;d1=abs<d>;whilefix<d1>>0d1=abs<d>/<2^m>;m=m+1;endbeq=fix<d1*2^b+.5>;beq=sign<d>.*beq.*2^<m-b-1>;第7章多采樣率信號處理例7-1在時域上顯示一個,信號頻率為0.042的正弦信號,然后以抽取因子3降采樣率,并在時域上顯示相應的結果,比較兩者在時域上的特點。解用MATLAB計算程序如下： M=3;%down-samplingfactor=3; fo=0.042;%signalfrequency=0.042; %generatetheinputsinusoidalsequence n=0:N-1; m=0:N*M-1; x=sin<2*pi*fo*