DSP課程設計報告-3

PAGEPAGE1課程設計(綜合實驗)報告名稱：數(shù)字信號處理課程設計院系：電子與通信工程系班級：學號：學生姓名：指導教師：設計周數(shù)：一周成績：日期：2010年12月20日

一、目的與要求掌握《數(shù)字信號處理基礎》課程的基本理論；掌握應用MATLAB進行數(shù)字信號處理的程序設計方法。二、主要內容設計題目及設計要求：已知低通數(shù)字濾波器的性能指標如下：，，，選擇合適的窗函數(shù)，設計滿足上述指標的數(shù)字線性相位FIR低通濾波器。用一個圖形窗口，包括四個子圖，分析顯示濾波器的單位沖激響應、相頻響應、幅頻響應和以dB為縱坐標的幅頻響應曲線。用雙線性變換法，設計滿足上述指標的數(shù)字ChebyshevI型低通濾波器。用一個圖形窗口，包括三個子圖，分析顯示濾波器的幅頻響應、以dB為縱坐標的幅頻響應和相頻響應。已知模擬信號其中，，，，取采樣頻率。以對進行取樣，得到。用一個圖形窗口，包括兩個子圖，分別顯示以及()的波形；用FFT對進行譜分析，要求頻率分辨率不超過5Hz。求出一個記錄長度中的最少點數(shù)，并用一個圖形窗口，包括兩個子圖，分別顯示以及的幅值；用要求1中設計的線性相位低通數(shù)字濾波器對進行濾波，求出濾波器的輸出，并用FFT對進行譜分析，要求頻率分辨率不超過5Hz。求出一個記錄長度中的最少點數(shù)，并用一個圖形窗口，包括四個子圖，分別顯示（）、、（）和的幅值；用要求2中設計的Chebyshev低通數(shù)字濾波器對進行濾波，求出濾波器的輸出，并用FFT對進行譜分析，要求頻率分辨率不超過5Hz。求出一個記錄長度中的最少點數(shù)，并用一個圖形窗口，包括四個子圖，分別顯示（）、、（）和的幅值。三、設計思路1）窗函數(shù)設計數(shù)字線性相位FIR低通濾波器1、根據(jù)設計給出的技術指標，計算允許的過渡帶寬度及阻帶衰減，初步選擇窗函數(shù)和N值。2、若選用理想低通逼近，則計算積分，求出理想低通的沖激響應hd(n);理想低通的截止頻率選擇為wc=(ws+wp)/2,hd(n)=sin[wc(n-a)]/pi*n-a);3、將hd(n)與窗函數(shù)相乘的FIR數(shù)字濾波器的沖激響應h(n)：h(n)=hd(n)w(n)4、計算FIR數(shù)字濾波器的頻率響應，并驗證是否達到所要求的指標。2）雙線性變換法設計的數(shù)字ChebyshevI型低通濾波器1、根據(jù)設計給出的數(shù)字濾波器的技術指標，計算模擬原型低通的指標2、模擬切比雪夫1型濾波器的計算，利用function[b,a]=afd_cheb1(wp,ws,Rp,As)進行切比雪夫1型模擬濾波器的設計；利用function[b,a]=cheb1ap_o(N,Rp,omega)非歸一化的切比雪夫1型模擬原型低通濾波器的設計3、雙線性變換法將Ha(s)轉換為H(z)4、直接型轉換為級聯(lián)型，計算數(shù)字濾波器的頻率響應，并驗證是否達到所要求的指標。3）模擬信號的波形產(chǎn)生和采樣及FFT頻譜分析和FIR濾波及IIR濾波N：一個記錄長度中的取樣數(shù)F：頻率分辨率fs：采樣周期最小記錄長度的選擇：由tp=NT=1/F可知，給定N=fs/F三、程序代碼1、主程序代碼%一、數(shù)字線性相位FIR低通濾波器clear;%數(shù)字濾波器的指標wp=0.26*pi;ws=0.41*pi;tr_width=ws-wp;%過渡帶寬N=ceil(6.6*pi/tr_width)+1;%按哈明窗計算所需的濾波器長度n=0:N-1;wc=(ws+wp)/2;%理想低通的截止頻率hd=ideal_lp(wc,N);%理想低通的沖激響應w_ham=(hamming(N))';%求窗函數(shù)h=hd.*w_ham;%實際沖激響應[db,mag,pha,grd,w]=freqz_m(h,[1]);%對設計結果進行檢驗delta_w=2*pi/1000;%頻率分辨率Rp=-(min(db(1:wp/delta_w+1)));%實際的通帶衰減As=-round(max(db(ws/delta_w+1:501)));%實際的最小阻帶衰減figure(1);%繪圖程序，分別顯示濾波器的單位沖激響應、相頻響應、幅頻響應和以dB為縱坐標的幅頻響應曲線subplot(2,2,1);stem(n,h);title('單位沖激響應');grid;axis([0N-1-0.10.5]);ylabel('沖激響應h(n)');xlabel('離散時間n');subplot(2,2,2);plot(w/pi,mag);title('幅頻響應');grid;axis([00.801.1]);ylabel('幅度');xlabel('以/pi為單位的頻率');subplot(2,2,3);plot(w/pi,db);title('以dB為縱坐標的幅頻響應');grid;axis([00.8-1000]);ylabel('對數(shù)幅度/dB');xlabel('以/pi為單位的頻率');subplot(2,2,4);plot(w/pi,pha);title('相頻響應');grid;axis([00.8-44]);ylabel('相位');xlabel('以/pi為單位的頻率');%二、雙線性變換法設計的數(shù)字ChebyshevI型低通濾波器Rp=0.75;As=50;%轉換為模擬原型低通指標T=1;Fs=1/T;omegap=(2/T)*tan(wp/2);omegas=(2/T)*tan(ws/2);%模擬切比雪夫1型濾波器的計算[cs,ds]=afd_cheb1(omegap,omegas,Rp,As);%雙線性變換法[b,a]=bilinear(cs,ds,Fs);%[C,B,A]=dir2cas(b,a);%直接型轉換為級聯(lián)型%計算數(shù)字濾波器的頻率響應[db,mag,pha,grd,w]=freqz_m(b,a);figure(2);subplot(2,2,1);plot(w/pi,mag);title('幅頻響應');grid;axis([00.801.1]);ylabel('幅度');xlabel('以/pi為單位的頻率');subplot(2,2,2);plot(w/pi,db);title('以dB為縱坐標的幅頻響應');grid;axis([01.0-30010]);ylabel('對數(shù)幅度/dB');xlabel('以/pi為單位的頻率');subplot(2,2,3);plot(w/pi,pha);title('相頻響應');grid;axis([00.8-44]);ylabel('相位');xlabel('以/pi為單位的頻率');%三（1）以fs=10khz取樣得到的x(n)(N=512)以及原始信號x(t)n=0:511;f1=120;f2=4980;f3=3250;f4=1150;fs=10000;T=1/fs;t=n*T;F=5;xn=2*sin(2*pi*f1*n*T)+5*sin(2*pi*f2*n*T)+8*cos(2*pi*f3*n*T)+7.5*cos(2*pi*f4*n*T);figure(3)subplot(2,1,1);plot(t,xn,'b');grid;axis([00.052-2525]);xlabel('t');ylabel('x(t)');title('連續(xù)信號x(t)');subplot(2,1,2);stem(n,xn,'b');grid;axis([0511-2525]);xlabel('n');ylabel('x(n)');title('取樣后的信號x(n)');%三（2）FFT對進行譜分析(N=2000)N=fs/F;n0=0:N-1;xn0=2*sin(2*pi*f1*n0*T)+5*sin(2*pi*f2*n0*T)+8*cos(2*pi*f3*n0*T)+7.5*cos(2*pi*f4*n0*T);XK=fft(xn0);m=[0:(length(xn0)-1)];figure(4)subplot(2,1,1);stem(n0,xn0,'b');grid;xlabel('n0');ylabel('x(n0)');title('取樣信號x(n0)');subplot(2,1,2);stem(m,abs(XK),'b');grid;xlabel('k');ylabel('X(k)');title('X(k)的幅值N=2000');%三（3）用要求1中設計的線性相位低通數(shù)字濾波器對x(n)進行濾波%y=filter(b,1,x);%可采用此函數(shù)給出濾波器的輸出y1=fftfilt(h,xn0);%y1=nonrec(h,xn0);%y1=conv(h,xn0);Ny1=[0:(length(y1)-1)];XK1=fft(xn0);AM1=abs(XK1);YK1=fft(y1);AM2=abs(YK1);figure(5)subplot(2,2,1);stem(n0,xn0);grid;xlabel('n0');ylabel('x(n0)');title('取樣信號x(n0)的幅值');subplot(2,2,3);stem(Ny1,y1);grid;xlabel('n0');ylabel('y1(n0)');title('濾波器的輸出y1(n0)的幅值');subplot(2,2,2);stem(n0,AM1);grid;xlabel('k');ylabel('X1(k)');title('X(k)的幅值');subplot(2,2,4);stem(Ny1,AM2);grid;xlabel('k');ylabel('Y1(k)');title('Y1(k)的幅值');%三（4）用要求2中設計的Chebyshev低通數(shù)字濾波器對進行濾波y2=filter(b,a,xn0);Ny2=[0:(length(y2)-1)];XK2=fft(xn0);AM1=abs(XK2);YK2=fft(y2);AM2=abs(YK2);figure(6)subplot(2,2,1);stem(n0,xn0);grid;xlabel('n0');ylabel('x(n0)');title('取樣信號x(n0)的幅值');subplot(2,2,3);stem(Ny2,y2);grid;xlabel('n0');ylabel('y2(n0)');title('濾波器的輸出y2(n0)的幅值');subplot(2,2,2);stem(n0,AM1);grid;xlabel('k');ylabel('X2(k)');title('X(k)的幅值');subplot(2,2,4);stem(Ny2,AM2);grid;xlabel('k');ylabel('Y2(k)');title('Y2(k)的幅值');2、所定義的函數(shù)代碼%hd=點0到N-1之間的理想沖激響應functionhd=ideal_lp(wc,M)%wc=截止頻率（弧度）%M=理想濾波器的長度(所選窗口的長度）alpha=(M-1)/2;n=0:M-1;m=n-alpha+eps;%加一個小數(shù)以避免零作除數(shù)hd=sin(wc*m)./(pi*m);%計算數(shù)字濾波器的頻率響應function[db,mag,pha,grd,w]=freqz_m(b,a);%db=[0到pi弧度]區(qū)間內的相對振幅(db)%mag=[0到pi弧度]區(qū)間內的絕對振幅%pha=[0到pi弧度]區(qū)間內的相位響應%grd=[0到pi弧度]區(qū)間內的群延遲%w=[0到pi弧度]區(qū)間內的501個頻率樣本向量%b=Ha(z)的分子多項式系數(shù)（對FIRb=h)%a=Ha(z)的分母多項式系數(shù)（對FIRa=[1])[H,w]=freqz(b,a,1000,'whole');H=(H(1:501))';w=(w(1:501))';mag=abs(H);%計算絕對振幅db=20*log10((mag+eps)/max(mag));%計算相對振幅(db)pha=angle(H);%計算相位響應grd=grpdelay(b,a,w);%計算群延遲%非歸一化的切比雪夫1型模擬原型低通濾波器的設計函數(shù)function[b,a]=cheb1ap_o(N,Rp,omega)[z,p,k]=cheb1ap(N,Rp);%設計一個階數(shù)為N，通帶紋波為Rp的歸一化的切比雪夫1型模擬原型低通濾波器的設計函數(shù)a=real(poly(p));aNn=a(N+1);p=p*omega;a=real(poly(p));aNu=a(N+1);k=k*aNu/aNn;b0=k;B=real(poly(z));b=k*B;%切比雪夫1型模擬濾波器的設計function[b,a]=afd_cheb1(wp,ws,Rp,As)ifwp<=0error('通帶必須大于0')endifws<=wperror('阻帶邊緣必須大于通帶邊緣')endif(Rp<=0)|(As<0)error('通帶波動或阻帶衰減必須大于0')endep=sqrt(10^(Rp/10)-1);A=10^(As/20);omegac=wp;omegar=ws/wp;g=sqrt(A*A-1)/ep;N=ceil(log10(g+sqrt(g*g-1))/log10(omegar+sqrt(o