DSP課程設(shè)計報告-3

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

一、目的與要求掌握《數(shù)字信號處理基礎(chǔ)》課程的基本理論；掌握應(yīng)用MATLAB進(jìn)行數(shù)字信號處理的程序設(shè)計方法。二、主要內(nèi)容設(shè)計題目及設(shè)計要求：已知低通數(shù)字濾波器的性能指標(biāo)如下：，，，選擇合適的窗函數(shù)，設(shè)計滿足上述指標(biāo)的數(shù)字線性相位FIR低通濾波器。用一個圖形窗口，包括四個子圖，分析顯示濾波器的單位沖激響應(yīng)、相頻響應(yīng)、幅頻響應(yīng)和以dB為縱坐標(biāo)的幅頻響應(yīng)曲線。用雙線性變換法，設(shè)計滿足上述指標(biāo)的數(shù)字ChebyshevI型低通濾波器。用一個圖形窗口，包括三個子圖，分析顯示濾波器的幅頻響應(yīng)、以dB為縱坐標(biāo)的幅頻響應(yīng)和相頻響應(yīng)。已知模擬信號其中，，，，取采樣頻率。以對進(jìn)行取樣，得到。用一個圖形窗口，包括兩個子圖，分別顯示以及()的波形；用FFT對進(jìn)行譜分析，要求頻率分辨率不超過5Hz。求出一個記錄長度中的最少點(diǎn)數(shù)，并用一個圖形窗口，包括兩個子圖，分別顯示以及的幅值；用要求1中設(shè)計的線性相位低通數(shù)字濾波器對進(jìn)行濾波，求出濾波器的輸出，并用FFT對進(jìn)行譜分析，要求頻率分辨率不超過5Hz。求出一個記錄長度中的最少點(diǎn)數(shù)，并用一個圖形窗口，包括四個子圖，分別顯示（）、、（）和的幅值；用要求2中設(shè)計的Chebyshev低通數(shù)字濾波器對進(jìn)行濾波，求出濾波器的輸出，并用FFT對進(jìn)行譜分析，要求頻率分辨率不超過5Hz。求出一個記錄長度中的最少點(diǎn)數(shù)，并用一個圖形窗口，包括四個子圖，分別顯示（）、、（）和的幅值。三、設(shè)計思路1）窗函數(shù)設(shè)計數(shù)字線性相位FIR低通濾波器1、根據(jù)設(shè)計給出的技術(shù)指標(biāo)，計算允許的過渡帶寬度及阻帶衰減，初步選擇窗函數(shù)和N值。2、若選用理想低通逼近，則計算積分，求出理想低通的沖激響應(yīng)hd(n);理想低通的截止頻率選擇為wc=(ws+wp)/2,hd(n)=sin[wc(n-a)]/pi*n-a);3、將hd(n)與窗函數(shù)相乘的FIR數(shù)字濾波器的沖激響應(yīng)h(n)：h(n)=hd(n)w(n)4、計算FIR數(shù)字濾波器的頻率響應(yīng)，并驗證是否達(dá)到所要求的指標(biāo)。2）雙線性變換法設(shè)計的數(shù)字ChebyshevI型低通濾波器1、根據(jù)設(shè)計給出的數(shù)字濾波器的技術(shù)指標(biāo)，計算模擬原型低通的指標(biāo)2、模擬切比雪夫1型濾波器的計算，利用function[b,a]=afd_cheb1(wp,ws,Rp,As)進(jìn)行切比雪夫1型模擬濾波器的設(shè)計；利用function[b,a]=cheb1ap_o(N,Rp,omega)非歸一化的切比雪夫1型模擬原型低通濾波器的設(shè)計3、雙線性變換法將Ha(s)轉(zhuǎn)換為H(z)4、直接型轉(zhuǎn)換為級聯(lián)型，計算數(shù)字濾波器的頻率響應(yīng)，并驗證是否達(dá)到所要求的指標(biāo)。3）模擬信號的波形產(chǎn)生和采樣及FFT頻譜分析和FIR濾波及IIR濾波N：一個記錄長度中的取樣數(shù)F：頻率分辨率fs：采樣周期最小記錄長度的選擇：由tp=NT=1/F可知，給定N=fs/F三、程序代碼1、主程序代碼%一、數(shù)字線性相位FIR低通濾波器clear;%數(shù)字濾波器的指標(biāo)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);%理想低通的沖激響應(yīng)w_ham=(hamming(N))';%求窗函數(shù)h=hd.*w_ham;%實際沖激響應(yīng)[db,mag,pha,grd,w]=freqz_m(h,[1]);%對設(shè)計結(jié)果進(jìn)行檢驗delta_w=2*pi/1000;%頻率分辨率Rp=-(min(db(1:wp/delta_w+1)));%實際的通帶衰減As=-round(max(db(ws/delta_w+1:501)));%實際的最小阻帶衰減figure(1);%繪圖程序，分別顯示濾波器的單位沖激響應(yīng)、相頻響應(yīng)、幅頻響應(yīng)和以dB為縱坐標(biāo)的幅頻響應(yīng)曲線subplot(2,2,1);stem(n,h);title('單位沖激響應(yīng)');grid;axis([0N-1-0.10.5]);ylabel('沖激響應(yīng)h(n)');xlabel('離散時間n');subplot(2,2,2);plot(w/pi,mag);title('幅頻響應(yīng)');grid;axis([00.801.1]);ylabel('幅度');xlabel('以/pi為單位的頻率');subplot(2,2,3);plot(w/pi,db);title('以dB為縱坐標(biāo)的幅頻響應(yīng)');grid;axis([00.8-1000]);ylabel('對數(shù)幅度/dB');xlabel('以/pi為單位的頻率');subplot(2,2,4);plot(w/pi,pha);title('相頻響應(yīng)');grid;axis([00.8-44]);ylabel('相位');xlabel('以/pi為單位的頻率');%二、雙線性變換法設(shè)計的數(shù)字ChebyshevI型低通濾波器Rp=0.75;As=50;%轉(zhuǎn)換為模擬原型低通指標(biāo)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);%直接型轉(zhuǎn)換為級聯(lián)型%計算數(shù)字濾波器的頻率響應(yīng)[db,mag,pha,grd,w]=freqz_m(b,a);figure(2);subplot(2,2,1);plot(w/pi,mag);title('幅頻響應(yīng)');grid;axis([00.801.1]);ylabel('幅度');xlabel('以/pi為單位的頻率');subplot(2,2,2);plot(w/pi,db);title('以dB為縱坐標(biāo)的幅頻響應(yīng)');grid;axis([01.0-30010]);ylabel('對數(shù)幅度/dB');xlabel('以/pi為單位的頻率');subplot(2,2,3);plot(w/pi,pha);title('相頻響應(yīng)');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對進(jìn)行譜分析(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è)計的線性相位低通數(shù)字濾波器對x(n)進(jì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中設(shè)計的Chebyshev低通數(shù)字濾波器對進(jìn)行濾波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=點(diǎn)0到N-1之間的理想沖激響應(yīng)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ù)字濾波器的頻率響應(yīng)function[db,mag,pha,grd,w]=freqz_m(b,a);%db=[0到pi弧度]區(qū)間內(nèi)的相對振幅(db)%mag=[0到pi弧度]區(qū)間內(nèi)的絕對振幅%pha=[0到pi弧度]區(qū)間內(nèi)的相位響應(yīng)%grd=[0到pi弧度]區(qū)間內(nèi)的群延遲%w=[0到pi弧度]區(qū)間內(nèi)的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);%計算相位響應(yīng)grd=grpdelay(b,a,w);%計算群延遲%非歸一化的切比雪夫1型模擬原型低通濾波器的設(shè)計函數(shù)function[b,a]=cheb1ap_o(N,Rp,omega)[z,p,k]=cheb1ap(N,Rp);%設(shè)計一個階數(shù)為N，通帶紋波為Rp的歸一化的切比雪夫1型模擬原型低通濾波器的設(shè)計函數(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型模擬濾波器的設(shè)計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