.16/161數字信號處理各章節重點知識的matlab實例第1,2章離散時間信號與系統例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給出了卷積結果的圖形,求得的結果存放在數組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個樣點的輸出,即該系統的單位脈沖響應。例1-3用MATLAB計算例1-2差分方程所對應的系統函數的DTFT。解例1-2差分方程所對應的系統函數為：其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對連續的單一頻率周期信號,按采樣頻率采樣,截取長度N分別選N=20和N=16,觀察其DFT結果的幅度譜。解此時離散序列,即k=8。用MATLAB計算并作圖,函數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結果,由于截取了兩個半周期,頻譜出現泄漏；<c>和<d>分別是N=16時的截取信號和DFT結果,由于截取了兩個整周期,得到單一譜線的頻譜。上述頻譜的誤差主要是由于時域中對信號的非整周期截斷產生的頻譜泄漏。例2-2用FFT計算兩個序列的互相關函數。解用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計算兩個序列的的互相關函數,其中x<n>={23521–100123530–1–2012}y<n>=x<n-4>+e<n>,e<n>為一隨機噪聲,在MATLAB中可以用隨機函數rand產生解用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>,我們看到最大值出現在m=4處,正好是y<n>對于x<n>的延遲。2.3<b>是x<n>的自相關函數,他和y<n>的區別除時間位置外,形狀也略不同,這是由于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設計一數字高通濾波器,它的通帶為400～500Hz,通帶內容許有0.5dB的波動,阻帶內衰減在小于317Hz的頻帶內至少為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計算的結果,可以看到模擬濾波器在Ω=∞處的三階零點通過高通變換后出現在ω=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計算的結果,可以看出數字濾波器將無窮遠點的二階零點映射為z=±1的二階零點,數字帶通濾波器的極點數是模擬低通濾波器的極點數的兩倍。例3-4一數字濾波器采樣頻率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數字濾波器,其指標為：通帶邊界頻率fc=800Hz,阻帶邊界fr=1000Hz,通帶波動阻帶最小衰減At=40dB,采樣頻率fs=4000Hz。解在MATLAB中可以用remezord和remez兩個函數設計,其結果如圖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'>函數remezord中的數組fedge為通帶和阻帶邊界頻率,數組mval是兩個邊界處的幅值,而數組dev是通帶和阻帶的波動,fs是采樣頻率單位為Hz。6章數字信號處理系統的實現例5-1求下列直接型系統函數的零、極點,并將它轉換成二階節形式解用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<'增益系數'>;disp<k>; sos=zp2sos<z,p,k>; disp<'二階節'>;disp<real<sos>>; zplane<num,den>輸入到"num"和"den"的分別為分子和分母多項式的系數。計算求得零、極點增益系數和二階節的系數：零點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增益系數1二階節1.0000-0.3885-0.55091.00001.15520.65111.00000.28850.36301.0000-1.05520.7679系統函數的二階節形式為：極點圖見圖5.1。例5-2分析五階橢圓低通濾波器的量化效應,其截止頻率為0.4,通帶紋波為0.4dB,最小的阻帶衰減為50dB。對濾波器進行截尾處理時,使用函數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表示系數是無限精度的理想濾波器的頻率響應〔以實線表示以及當濾波器系數截尾到5位時的頻率響應〔以短線表示。由圖可知,系數量化對頻帶的邊緣影響較大,經系數量化后,增加了通帶的波紋幅度,減小了過渡帶寬,并且減小了最小的阻帶衰減。圖5.1〔b給出了系數量化以前和系數量化以后的橢圓低通濾波器的零極點位置。由圖可知,系數的量化會使零極點的位置與它們的理想的標稱位置相比發生顯著的改變。在這個例子中,靠近虛軸的零點的位置變動最大,并且移向靠它最近的極點的位置。只要對程序稍作改變就可以分析舍入量化的影響。為了研究二進制數量化效應對數字濾波器的影響,首先需要將十進制表示的濾波器系數轉換成二進制數并進行量化,二進制數的量化既可以通過截尾法也可以通過舍入法實現。我們提供了如下的兩個MATLAB程序：a2dT.m和a2dR.m,這兩段程序分別將向量d中的每一個數按二進制數進行截尾或舍入量化,量化的精度是小數點以后保留b位,量化后返回的向量為beq。functionbeq=a2dT<d,b>%beq=a2dT<d,b>將十進制數利用截尾法得到b位的二進制數,%然后將該二進制數再轉換為十進制數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>將十進制數利用舍入法得到b位的二進制數%然后將該二進制數再轉換為十進制數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*