華北電力大學(xué)科技學(xué)院數(shù)字信號(hào)處理課程設(shè)計(jì)

1、科 技 學(xué) 院課程設(shè)計(jì)(綜合實(shí)驗(yàn))報(bào)告( 2015- 2016年度第一學(xué)期)名 稱：數(shù)字信號(hào)處理課程設(shè)計(jì) 題 目：MATLAB編程 院 系：信息工程系 班 級(jí)： 13K2 學(xué) 號(hào)：31 學(xué)生姓名： 指導(dǎo)教師：孫老師 設(shè)計(jì)周數(shù)： 2 成 績(jī)： 日期：2015年 12月 18日一、課程設(shè)計(jì)(綜合實(shí)驗(yàn))的目的與要求一、 目的與要求1 掌握數(shù)字信號(hào)處理基礎(chǔ)課程的基本理論；2 掌握應(yīng)用MATLAB進(jìn)行數(shù)字信號(hào)處理的程序設(shè)計(jì)方法。二、 主要內(nèi)容設(shè)計(jì)題目及設(shè)計(jì)要求：已知低通數(shù)字濾波器的性能指標(biāo)如下：，要求：1. 選擇合適的窗函數(shù)，設(shè)計(jì)滿足上述指標(biāo)的數(shù)字線性相位FIR低通濾波器。用一個(gè)圖形窗口，包括四個(gè)子圖，

2、分析顯示濾波器的單位沖激響應(yīng)、相頻響應(yīng)、幅頻響應(yīng)和以dB為縱坐標(biāo)的幅頻響應(yīng)曲線。2. 用雙線性變換法，設(shè)計(jì)滿足上述指標(biāo)的數(shù)字Chebyshev I型低通濾波器。用一個(gè)圖形窗口，包括三個(gè)子圖，分析顯示濾波器的幅頻響應(yīng)、以dB為縱坐標(biāo)的幅頻響應(yīng)和相頻響應(yīng)。3. 已知模擬信號(hào) 其中，取采樣頻率。要求：（1）以對(duì)進(jìn)行取樣，得到。用一個(gè)圖形窗口，包括兩個(gè)子圖，分別顯示以及()的波形；（2）用FFT對(duì)進(jìn)行譜分析，要求頻率分辨率不超過(guò)5Hz。求出一個(gè)記錄長(zhǎng)度中的最少點(diǎn)數(shù)，并用一個(gè)圖形窗口，包括兩個(gè)子圖，分別顯示以及的幅值；（3）用要求1中設(shè)計(jì)的線性相位低通數(shù)字濾波器對(duì)進(jìn)行濾波，求出濾波器的輸出，并用FFT對(duì)

3、進(jìn)行譜分析，要求頻率分辨率不超過(guò)5Hz。求出一個(gè)記錄長(zhǎng)度中的最少點(diǎn)數(shù)，并用一個(gè)圖形窗口，包括四個(gè)子圖，分別顯示（）、（）和的幅值；（4）用要求2中設(shè)計(jì)的Chebyshev低通數(shù)字濾波器對(duì)進(jìn)行濾波，求出濾波器的輸出，并用FFT對(duì)進(jìn)行譜分析，要求頻率分辨率不超過(guò)5Hz。求出一個(gè)記錄長(zhǎng)度中的最少點(diǎn)數(shù)，并用一個(gè)圖形窗口，包括四個(gè)子圖，分別顯示（）、（）和的幅值。備注：(1) 要求編寫(xiě)一個(gè)主程序，完成上述問(wèn)題。(2) 要求各結(jié)果圖均標(biāo)出圖題以及橫縱軸的名稱。(3) 要求給主要語(yǔ)句加上標(biāo)注。二、設(shè)計(jì)（實(shí)驗(yàn)）正文1：設(shè)計(jì)思路:1）窗函數(shù)設(shè)計(jì)數(shù)字線性相位FIR低通濾波器1、根據(jù)設(shè)計(jì)給出的技術(shù)指標(biāo)，計(jì)算允許的過(guò)

4、渡帶寬度及阻帶衰減，初步選擇窗函數(shù)和N值(可選用hamming窗，Blackman窗和Kaiser窗,這里選用Hamming窗)。2、若選用理想低通逼近，則計(jì)算積分，求出理想低通的沖激響應(yīng)hd(n);理想低通的截止頻率選擇為wc=(ws+wp)/2,hd(n)=sinwc(n-a)/pi*n-a);3、將hd(n)與窗函數(shù)相乘的FIR數(shù)字濾波器的沖激響應(yīng)h(n)：h(n)=hd(n)w(n)4、計(jì)算FIR數(shù)字濾波器的頻率響應(yīng)，并驗(yàn)證是否達(dá)到所要求的指標(biāo)。2）雙線性變換法設(shè)計(jì)的數(shù)字Chebyshev I型低通濾波器1、根據(jù)設(shè)計(jì)給出的數(shù)字濾波器的技術(shù)指標(biāo)，計(jì)算模擬原型低通的指標(biāo)2、模擬切比雪夫1型

5、濾波器的計(jì)算，利用functionb,a=afd_cheb1(wp,ws,Rp,As)進(jìn)行切比雪夫1型模擬濾波器的設(shè)計(jì)；利用functionb,a=cheb1ap_o(N,Rp,omega)進(jìn)行非歸一化的切比雪夫1型模擬原型低通濾波器的設(shè)計(jì)3、雙線性變換法將Ha(s)轉(zhuǎn)換為H(z)4、直接型轉(zhuǎn)換為級(jí)聯(lián)型，計(jì)算數(shù)字濾波器的頻率響應(yīng)，并驗(yàn)證是否達(dá)到所要求的指標(biāo)。3）模擬信號(hào)的波形產(chǎn)生和采樣及FFT頻譜分析和FIR濾波及IIR濾波N：一個(gè)記錄長(zhǎng)度中的取樣數(shù)F：頻率分辨率 fs：采樣周期最小記錄長(zhǎng)度的選擇：由 tp=NT=1/F可知，給定抽樣頻率和頻率分辨率后，就可以得出最小記錄長(zhǎng)度 N=fs/Fff

6、t(x)：計(jì)算N點(diǎn)的DFT,N是序列的長(zhǎng)度f(wàn)ft(x，L)：計(jì)算L點(diǎn)的DFT,若N>L，則將原序列截短為L(zhǎng)點(diǎn)的序列，再計(jì)算其L點(diǎn)的DFT；若N<L，則將原序列補(bǔ)零為L(zhǎng)點(diǎn)的序列，再計(jì)算其L點(diǎn)的DFT；ifft(x)：計(jì)算N點(diǎn)的IDFT,N是頻域序列的長(zhǎng)度ifft(x，L)：計(jì)算L點(diǎn)的IDFT,若N>L，則將原序列截短為L(zhǎng)點(diǎn)的序列，再計(jì)算其L點(diǎn)的IDFT；若N<L，則將原序列補(bǔ)零為L(zhǎng)點(diǎn)的序列，再計(jì)算其L點(diǎn)的IDFT； 為了提高fft和ifft的計(jì)算效率，應(yīng)盡量使序列長(zhǎng)度N為2的冪次方，或?qū)⑿蛄醒a(bǔ)零使其長(zhǎng)度為2的冪次方。filter(b,a,x)：可采用此函數(shù)給出濾波器的輸出

7、b=H(z)的分子多項(xiàng)式系數(shù)（對(duì)FIR b=h)a=H(z)的分母多項(xiàng)式系數(shù)（對(duì)FIR a=1)2：主函數(shù)代碼以及需調(diào)用的函數(shù)（1）主程序代碼(以下代碼中，由于復(fù)制過(guò)來(lái)，word不認(rèn)識(shí)MATLAB中的文字，所以自己又重新打了，才會(huì)有的是沒(méi)顏色，但并沒(méi)有影響)%1wp=0.26*pi;ws=0.41*pi;As=50;Rp=0.75;tr_width=ws-wp;N=81%未采用經(jīng)驗(yàn)公式，若需要用的同學(xué)把下面的注釋去掉;用80也可以%N=ceil(6.6*pi/tr_width)+1%經(jīng)驗(yàn)公式,求采樣點(diǎn)數(shù)n=0:1:N-1;wc=(ws+wp)/2; %理想低通的截止頻率hd=ideal_lp(

8、wc,N); %理想低通的沖激響應(yīng)w_ham=(hamming(N)' %哈明窗的窗函數(shù)h=hd.*w_ham; %FIR濾波器的實(shí)際沖激響應(yīng)db,mag,pha,grd,w=freqz_m(h,1); %計(jì)算數(shù)字濾波器的頻率響應(yīng)delta_w=2*pi/1000; %取1000點(diǎn)Rp=-(min(db(1:1:wp/delta_w+1) %實(shí)際的通帶衰減,間隔1dbAs=-round(max(db(ws/delta_w+1:1:501) %實(shí)際的最小阻帶衰減figure(1); % 分別顯示濾波器的單位沖激響應(yīng)、相頻響應(yīng)、幅頻響應(yīng)和以dB為縱坐標(biāo)的幅頻響應(yīng)曲線。subplot(2,2

9、,1);stem(n,hd,'g');title('(理想)單位沖激響應(yīng)');grid;ylabel('hd(n)');xlabel('n');subplot(2,2,2);plot(w/pi,pha,'r');title('相位響應(yīng)');grid;ylabel('相位');xlabel('以pi為單位的頻率');subplot(2,2,3);plot(w/pi,mag,'k');title('幅度響應(yīng)');grid;ylabel(&#

10、39;|幅度|');xlabel('以pi為單位的頻率');subplot(2,2,4);plot(w/pi,db,'b');title('幅度響應(yīng)(dB)');grid;ylabel('對(duì)數(shù)幅度dB');xlabel('以pi為單位的頻率');可以對(duì)其進(jìn)行優(yōu)化%2wp=0.26*pi;ws=0.41*pi;Rp=0.75;As=50;%轉(zhuǎn)換為模擬濾波器的指標(biāo)（模擬原型低通指標(biāo)）T=1;Fs=1/T;omegap=(2/T)*tan(wp/2);omegas=(2/T)*tan(ws/2);cs,ds=af

11、d_cheb1(omegap,omegas,Rp,As); %模擬切比雪夫I型低通濾波器的計(jì)算,求解拉普拉斯變換的分子與分母多項(xiàng)式的系數(shù)數(shù)組b,a=bilinear(cs,ds,Fs); %雙線性變換法db,mag,pha,grd,w=freqz_m(b,a); %計(jì)算數(shù)字濾波器的頻率響應(yīng),頻域上采樣并記錄幅頻和相頻響應(yīng)%分別顯示數(shù)字切比雪夫I型低通濾波器的幅頻響應(yīng)、以dB為縱坐標(biāo)的幅頻響應(yīng)和相頻響應(yīng)figure(2)subplot(3,1,1);plot(w/pi,mag,'r');title('幅度響應(yīng)');grid;ylabel('|幅度|'

12、;);xlabel('以pi為單位的頻率');subplot(3,1,2);plot(w/pi,db,'g');title('以dB為單位的幅度響應(yīng)');grid;ylabel('對(duì)數(shù)幅度/dB');xlabel('以pi為單位的頻率');subplot(3,1,3);plot(w/pi,pha,'b');title('相位響應(yīng)');grid;ylabel('相位');xlabel('以pi為單位的頻率');可以對(duì)其進(jìn)行優(yōu)化%3N=0:511; %做5

13、12點(diǎn)DFTfs=10000; %采樣頻率t=0:0.0001:1;xt=2*sin(2*pi*120*t)+5*sin(2*pi*4980*t)+8*cos(2*pi*3250*t)+7.5*cos(2*pi*1150*t);T=1/fs; %采樣周期t0=N*T; %時(shí)間離散化Xn=2*sin(2*120*pi*t0)+5*sin(2*4980*pi*t0)+8*cos(2*pi*3250*t0)+7.5*cos(2*pi*1150*t0);figure(3) %顯示x(n)(0<=n<=511)以及x(t)的波形；subplot(2,1,1);stem(N,Xn,'r

14、');title('Xn');grid;ylabel('Xn');xlabel('n');subplot(2,1,2);plot(t,xt,'b');title('xt曲線');grid;ylabel('xt');xlabel('t');可以對(duì)其進(jìn)行優(yōu)化%4F=5; %頻率分辨率Nx=fs/F; %取樣數(shù)(因?yàn)榇藭r(shí)fs遠(yuǎn)大于x(t)中的最高頻率,應(yīng)該用fs才能分辨出四個(gè)譜峰f1f4),計(jì)算得Nx=2000點(diǎn)T=1/fs; %采樣周期%n=0:Nx-1; %做2000點(diǎn)DFTt

15、0=n*T;Xn=2*sin(2*120*pi*t0)+5*sin(2*4980*pi*t0)+8*cos(2*pi*3250*t0)+7.5*cos(2*pi*1150*t0);k=0:Nx-1; Xk=fft(Xn); AmX=abs(Xk); figure(4) % 顯示x(n)以及X(k)的幅值；subplot(2,1,1);stem(n,Xn,'g');title('Xn');grid;ylabel('Xn');xlabel('n');subplot(2,1,2);stem(k,AmX,'k');titl

16、e('Xk');grid;ylabel('AmX');xlabel('k');可以對(duì)其進(jìn)行優(yōu)化%5Fy1=5; %頻率分辨率Nx=fs/Fy1;Xn=2*sin(240*pi*t0)+5*sin(2*4980*pi*t0)+8*cos(2*pi*3250*t0)+7.5*cos(2*pi*1150*t0);k=0:Nx-1; Xk=fft(Xn); magX=abs(Xk);y1n=filter(h,1,Xn); %利用遞歸濾波器(IIR)或非遞歸濾波器(FIR)對(duì)數(shù)據(jù)進(jìn)行濾波(對(duì)FIR b=h,a=1),Xn經(jīng)過(guò)用哈明窗設(shè)計(jì)的線性相位低通濾波器

17、濾波得到y(tǒng)1nYk1=fft(y1n); %對(duì)y1n進(jìn)行FFT變換AmYk1=abs(Yk1); figure(5) %顯示x(n) (0<=n<=Nx-1)、X(k)、y1(n) (0<=n<=Ny1-1)和Y1(k)的幅值；subplot(2,2,1);stem(n,Xn,'b');title('Xn');grid;ylabel('Xn');xlabel('n');subplot(2,2,2);stem(k,magX,'r');title('Xk');grid;ylabe

18、l('magX');xlabel('k');subplot(2,2,3);stem(n,y1n,'g');title('y1n');grid;ylabel('y1n');xlabel('n');subplot(2,2,4);stem(k,AmYk1,'y');title('Yk1');grid;ylabel('magY1');xlabel('k');可以對(duì)其進(jìn)行優(yōu)化%6Fy2=5; %頻率分辨率Nx=fs/Fy1;Xn=2*sin(2*

19、120*pi*t0)+5*sin(2*4980*pi*t0)+8*cos(2*pi*3250*t0)+7.5*cos(2*pi*1150*t0);y2n=filter(b,a,Xn); %利用遞歸濾波器(IIR)或非遞歸濾波器(FIR)對(duì)數(shù)據(jù)進(jìn)行濾波(b和a值在要求2中已求出，直接調(diào)用)XK2=fft(Xn); %對(duì)Xn進(jìn)行DFT運(yùn)算AM1=abs(XK2); %取模值|XK| YK2=fft(y2n); %對(duì)y2n進(jìn)行DFT運(yùn)算AM2=abs(YK2); %取模值|YK|figure(6) %顯示x(n)(0<=n<=Nx-1)、X(k)、y2(n)(0<=n<=Ny

20、2-1)和Y2(k)的幅值subplot(2,2,1);stem(n,Xn,'b');grid;xlabel('n');ylabel('x(n)');title('取樣信號(hào)x(n)');subplot(2,2,2);stem(n,AM1,'r');grid;xlabel('k');ylabel('X(k)');title('X(k)µÄ·ùÖµ');subplot(2,2,3);stem(n,y2n,'

21、;g');grid;xlabel('n');ylabel('y2(n)');title('濾波器的輸出y2(n)');subplot(2,2,4);stem(n,AM2,'k');grid;xlabel('k');ylabel('Y2(k)');title(' Y2(k)的幅值');可以對(duì)其進(jìn)行優(yōu)化(2)需要調(diào)用的函數(shù)代碼1)functionb,a=afd_cheb1(wp,ws,Rp,As)%if wp<=0 error('通帶必須大于0')endif

22、 ws<=wp error('阻帶邊緣必須大于通帶邊緣 ')endif(Rp<=0)|(As<0) error('通帶波動(dòng)或阻帶衰減必須大于0 ')end%ep=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(omegar*omegar-1);fprintf('n*切比雪夫?yàn)V波器的階數(shù)=%2.0fn',N);b,a=cheb1ap_o(N,Rp,om

23、egac);2）function b,a=bilinear(cs,ds,Fs)R,p,k=residue(cs,ds);p=exp(p*Fs);b,a=residuez(R,p,k);b=real(b');a=real(a');3）functionb,a=cheb1ap_o(N,Rp,omega)z,p,k=cheb1ap(N,Rp);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;4）functionb0,B,A=dir2cas(b,a)b0=b(1);b=b/b0;a0=a(1);a=a/a0;b0=b0/a0;%M=length(b);N=length(a);if N>M b=b zeros(1,N-M); elsif M>N a=a zeros(1,M-N); N=M;else NM=0;end%K=floor(N/2);B=zeros(K,3);A=zeros(K,3);if K*2=N; b=b 0; a=