數(shù)字信號處理實驗報告

數(shù)字信號處理實驗報告專業(yè)：姓名：學號：指導教師：實驗一序列、頻譜、卷積實驗目的掌握序列的輸入方法；熟悉不同序列的特征；了解確定性信號譜分析的方法；驗證卷積的計算過程；實驗要求利用matlab程序，生成幾種常用的序列，如矩形序列，單位脈沖序列；繪制圖形，觀察序列特征；研究其頻率特性，繪制圖形，觀察頻率響應特征；利用matlab程序，驗證卷積的過程；實驗步驟矩形序列（1）生成長度為N的矩形序列，觀察并記錄生成的圖形；n=1:50x=sign(sign(10-n)+1);closeall;subplot(3,1,1);stem(x);title('單位矩形信號序列');（2）研究其頻率特性，，分別研究其幅頻特性和相頻特性，觀察并記錄生成的圖形；k=-25:25;X=x*(exp(-j*pi/25)).^(n'*k);magX=abs(X);%繪制x(n)的幅度譜subplot(3,1,2);stem(magX);title('單位矩形信號的幅度譜');angX=angle(X);%繪制x(n)的相位譜subplot(3,1,3);stem(angX);title('單位矩形信號的相位譜')答：根據(jù)以上兩個源程序得到的圖形如下圖1.1：圖1.1分析：幅度譜可以發(fā)現(xiàn)其幅度呈偶對稱，而相位譜與幅度普均為離散信號。其頻率特性就是通過傅里葉變換得來的，是用頻率信息表示時域信息，與連續(xù)信號的相同的矩形函數(shù)的頻率特性對應，只不過一個是離散的一個是連續(xù)的。單位脈沖序列生成單位脈沖序列，觀察并記錄生成的圖形；n=1:50;%定義序列的長度是50x=zeros(1,50);%注意：MATLAB中數(shù)組下標從1開始x(1)=1;closeall;subplot(3,1,1);stem(x);title('單位沖擊信號序列');研究其頻率特性，，分別研究其幅頻特性和相頻特性，觀察并記錄生成的圖形；k=-25:25;X=x*(exp(-j*pi/12.5)).^(n'*k);magX=abs(X);%繪制x(n)的幅度譜subplot(3,1,2);stem(magX);title('單位沖擊信號的幅度譜');angX=angle(X);%繪制x(n)的相位譜subplot(3,1,3);stem(angX);title('單位沖擊信號的相位譜')答：根據(jù)以上兩個源程序得到的圖形如下圖1.2：圖1.2分析：由上圖觀察單位沖激信號的頻率響應可知，其幅度譜為矩形函數(shù)，相位譜呈奇對稱。與連續(xù)的沖激響應對應，只是這里是有限序列。卷積過程，n=1:50;%定義序列的長度是50hb=zeros(1,50);%注意：MATLAB中數(shù)組下標從1開始hb(1)=1;hb(2)=2.5;hb(3)=2.5;hb(4)=1;closeall;subplot(3,1,1);stem(hb);title('系統(tǒng)hb[n]');m=1:50;%定義序列的長度是50A=444.128;%設置信號有關的參數(shù)a=50*sqrt(2.0)*pi;T=0.001;%采樣率w0=50*sqrt(2.0)*pi;x=A*exp(-a*m*T).*sin(w0*m*T);subplot(3,1,2);stem(x);title('輸入信號x[n]');y=conv(x,hb);subplot(3,1,3);stem(y);title('輸出信號y[n]');答：根據(jù)以上兩個源程序得到的圖形如下圖1.3：圖1.3分析：離散信號的卷積為相對移動然后累計相加，所以卷積的結(jié)果最高可達到1000.思考：如何生成實指數(shù)序列？寫出對應的matlab程序；答：（1）程序如下所示：a1=2n=1:50x1=(a1.^n)subplot(1,1,1)stem(x1);title('實指數(shù)序列')所得圖形如下圖1.4所示：編寫程序驗證卷積定律。答：編寫的源程序如下：x=[12031];%原序列y=[12031];%直接得到結(jié)果z=conv(x,y);figure(1),subplot(311),stem(x),title('直接法x的序列');axis([1904]);subplot(312),stem(y),title('直接法y的序列');axis([1904]);subplot(313),stem(z),title('直接法得到的卷積z的序列');axis([19020]);%FFTN=10;x1=[xzeros(1,N-length(x))];y1=[yzeros(1,N-length(y))];X1=fft(x1);Y1=fft(y1);Z1=X1.*Y1;z1=ifft(Z1);figure(2),subplot(321),stem(x1),title('x序列');subplot(322),stem(real(X1)),title('FFT法的傅里葉變換實部x');subplot(323),stem(y1),title('y序列');subplot(324),stem(real(Y1)),title('FFT法的傅里葉變換實部y');subplot(325),stem(z1),title('z序列');subplot(326),stem(real(Z1)),title('傅里葉變換后相乘結(jié)果z');N=6;x2=[xzeros(1,N-length(x))];y2=[yzeros(1,N-length(y))];X2=fft(x2);Y2=fft(y2);Z2=X2.*Y2;z2=ifft(Z2);figure(3),subplot(321),stem(x1),title('x序列2');subplot(322),stem(real(X1)),title('FFT法的傅里葉變換實部x2');subplot(323),stem(y1),title('y序列2');subplot(324),stem(real(Y1)),title('FFT法的傅里葉變換實部y2');subplot(325),stem(z1),title('z序列');subplot(326),stem(real(Z1)),title('傅里葉變換后相乘結(jié)果z2');得到驗證如下圖1.4至1.6：圖1.4圖1.5圖1.6分析：觀察以上三個圖可以發(fā)現(xiàn)驗證了卷積定理，即時域的卷積等于頻域的相乘算式為f(t)*h(t)=F(jw)H(jw)。與直接用MATLAB內(nèi)部卷積函數(shù)得到的結(jié)果也一樣，并且對比圖1.5和1.6可知快速傅里葉變換所用的點數(shù)的多少不影響卷積的結(jié)果。實驗二基2－FFT算法編程實驗目的掌握離散傅立葉變換的原理；熟悉基2－FFT算法流程，學會編寫算法程序；認識到基2－FFT算法對離散傅立葉運算的速度提高的重要性；實驗要求自己編寫FFT算法程序；調(diào)試算法程序；測試采用FFT算法后，運算量和運算時間的減少量；實驗步驟編寫matlab程序，完成輸入倒位序，輸出自然順序，按照時間抽選的基2－FFT算法程序的編寫；答：算法源程序如下：clc;closeall;formatcompact;xn=[1,2,3,4,5,6,7];%可取任意序列M=nextpow2(length(xn));N=2^M;form=0:N/2-1;%旋轉(zhuǎn)因子指數(shù)范圍WN(m+1)=exp(-j*2*pi/N)^m;%計算旋轉(zhuǎn)因子endA=[xn,zeros(1,N-length(xn))];%數(shù)據(jù)輸入disp('輸入到各存儲單元的數(shù)據(jù)：'),disp(A)%數(shù)據(jù)倒序操作J=0;%給倒序數(shù)賦初值forI=0:N-1;%按序交換數(shù)據(jù)和算倒序數(shù)ifI<J;%條件判斷及數(shù)據(jù)交換T=A(I+1);A(I+1)=A(J+1);A(J+1)=T;end%算下一個倒序數(shù)K=N/2;whileJ>=K;J=J-K;K=K/2;endJ=J+K;enddisp('倒序后各存儲單元的數(shù)據(jù)'),disp(A);%分級按序一次進行蝶形運算forL=1:M;%分級計算disp('運算級次'),disp(L);B=2^(L-1);forR=0:B-1;%各級按序碟算P=2^(M-L)*R;forK=R:2^L:N-2;%每序依次計算T=A(K+1)+A(K+B+1)*WN(P+1);A(K+B+1)=A(K+1)-A(K+B+1)*WN(P+1);A(K+1)=T;endenddisp('本級運算后各存儲單元的數(shù)據(jù)'),disp(A);enddisp('輸出各存儲單元的數(shù)據(jù)'),Xk=A,disp('調(diào)用FFT函數(shù)運算的結(jié)果'),fftxn=fft(xn,N),調(diào)試程序；答：調(diào)試上題程序得到結(jié)果如下圖2.1至2.2：圖2.1圖2.2已知，利用編寫的程序求；答：將源程序第二行改為：xn=0:1:31;運行可得結(jié)果如下圖2.3至2.4：圖2.3圖2.4答：源程序如下：clc;closeall;n=[0:31];x=n+1;h=[2,-1,1];y=hsolpsav(x,h,N)其中hsolpsav函數(shù)定義如下：function[y]=hsolpsav(x,h,N)%用fft的重疊保留法做分段卷積%%[y]=hsolpsav(x,h,N)%y=輸出序列%x=輸入序列%h=單位沖擊響應%N=分段長必須是二的冪次%N=2^(ceil(log10(N)/log10(2)));Lx=length(x);M=length(h);M1=M-1;L=N-M1;H=fft(h,N);%K=floor(((Lx+M1)/L)+1);%分段數(shù)目x=[zeros(1,M1),x,zeros(1,N-1)];%前置M1個零點Y=zeros(K,N);fork=0:K-1xk=fft(x(k*L+1:k*L+N));Y(k+1,:)=real(ifft(xk.*H));endY=Y(:,M:N)';y=(Y(:))';%去掉前面M1個值運行結(jié)果如下圖2.5：圖2.5答：源程序如下：clc;clear;closeall;N=input('N=');T=0.025;n=1:N;%原始數(shù)據(jù)D=2*pi/(N*T);%計算分辨率xa=cos(20*n*T);%有限長余弦序列Xa=T*fftshift(fft(xa,N));Xa(1)%求x（n）的DFT，并移到對稱位置k=floor(-(N-1)/2:(N-1)/2);%對w=0對稱的奈奎斯特頻率下標向量TITLE=sprintf('N=%i,L=%i,N,N*T');plot(k*D,abs(Xa));axis([-30,30,0,max(abs(Xa))+2]);xlabel('omega');ylabel('|X(j\Omega)|')title(TITLE)運行結(jié)果如下圖2.6至2.8：圖2.6N=100時得到的頻譜圖圖2.7N=500時得到的頻譜圖圖2.8N=1000時得到的頻譜圖分析：由以上三個圖觀察分析可知：當采樣點越多，得到的頻譜圖越精確越接近連續(xù)信號的理想結(jié)果。思考如何編寫重疊相加法算法程序？答：源程序如下（主程序）：clc;closeall;n=[0:31];x=n+1;h=[2,-1,1];y=overlap_add(x,h,4)其中overlap_add函數(shù)定義如下：%%重疊相加法，重疊保留法實現(xiàn)長短序列的線性卷積functiony=overlap_add(x,h,M)%x為較長的輸入序列，h為較短的系統(tǒng)函數(shù)，M為分段大小%前期處理N=length(h);%序列h(n)的長度Lx=length(x);%序列x(n)的長度ifN>M%算法要求N<=MM=N+1;endL=M+N-1;%用每段循環(huán)卷積計算線性卷積所需點數(shù)T=ceil(Lx/M);%分段數(shù)，ceil向上取整t=zeros(1,N-1);%緩存序列初始化x=[x,zeros(1,(T+1)*M-Lx)];%最后一個不足M的分段補零y=zeros(1,(T+1)*M);%生成輸出序列y(n)%核心算法fori=0:Txi=i*M+1;x_seg=x(xi:xi+M-1);%低點數(shù)計算時的分段x(n)y_seg=lin_conv(x_seg,h,L);%計算i分段和h的循環(huán)線性卷積y_seg(1:N-1)=y_seg(1:N-1)+t(1:N-1);%與前一段卷積的后N-1位重疊相加t(1:N-1)=y_seg(M+1:L);%緩存序列更新y(xi:xi+M-1)=y_seg(1:M);%每卷積一段輸出M個點endy=y(1:Lx+N-1);%最終輸出序列End其中overlap_add函數(shù)調(diào)用的函數(shù)lin_conv定義如下：%%利用FFT循環(huán)卷積計算線性卷積的函數(shù)functiony=lin_conv(x1,x2,L)%利用FFT循環(huán)卷積計算線性卷積y1=fft(x1,L);y2=fft(x2,L);Yk=y1.*y2;y=ifft(Yk);end運行結(jié)果如下圖2.9：圖2.9分析：對比用函數(shù)hsolpsav函數(shù)相比結(jié)果，重疊相加法可以減少運算次數(shù)，使卷積速度很快，而且沒有不必要的0結(jié)果。編寫完成的基2－FFT算法程序，如何變成基2－IFFT算法程序？答：報告老師，由于能力有限還有時間緊迫再加上電腦太卡了，就沒有做某些思考題，希望老師諒解！實驗三IIR數(shù)字濾波器的設計一、實驗目的掌握IIR數(shù)字濾波器設計的方法；熟悉軟件編程實現(xiàn)的方法；認識到數(shù)字濾波器在數(shù)字信號處理理論中的重要性；了解IIR數(shù)字濾波器實現(xiàn)方法；二、實驗要求了解和掌握各種數(shù)字濾波器的設計方法；整理低通數(shù)字濾波器的設計算法和設計步驟；編寫和調(diào)試程序；三、實驗步驟編寫程序，以巴特沃思模擬低通濾波器為樣本濾波器，設計基于雙線性變換法的低通數(shù)字低通濾波器；答：源程序如下：%數(shù)字濾波器指標：wp=0.2*pi;%數(shù)字通帶頻率（Hz）ws=0.3*pi;%數(shù)字阻帶頻率（Hz）Rp=1;%通帶波動（dB）As=15;%阻帶波動（dB）%模擬原型指標的頻率逆映射Fs=10000;T=1/Fs;%抽樣頻率Fs（Hz）OmegaP=(2/T)*tan(wp/2);%原型通帶頻率預修正OmegaS=(2/T)*tan(ws/2);%原型阻帶頻率預修正ep=sqrt(10^(Rp/10)-1);%通帶波動參數(shù)Ripple=sqrt(1/(1+ep*ep));%通帶波Attn=1/(10^(As/20));%阻帶衰減%模擬巴特沃思原型濾波器計算：[N,OmegaC]=buttord(OmegaP,OmegaS,Rp,As,'s');%原型的階數(shù)和截止頻率的計算[z0,p0,k0]=buttap(N);%歸一化巴特沃思原型設計函數(shù)p=p0*OmegaC;z=z0*OmegaC;%將零極點乘以OmegaC，得到非歸一化零極點k=k0*OmegaC^N;%將k0乘以OmegaC^N，得到非歸一化kba0=real(poly(z0));ba0=k0*ba0;%由零點計算分子系數(shù)向量aa0=real(poly(p0))%由極點計算分母系數(shù)向量ba=real(poly(z));ba=k*ba;%由零點計算分子系數(shù)向量aa=real(poly(p))%由極點計算分母系數(shù)向量[bd,ad]=bilinear(ba,aa,Fs);%調(diào)用雙線性變換法函數(shù)[bd1,ad1]=bilinear(ba0,aa0,Fs/OmegaC);%雙線性變換法%繪圖figure(1);subplot(1,1,1)[H,w]=freqz(bd,ad,1000,'whole');%計算數(shù)字系統(tǒng)頻率響應H=(H(1:1:501))';w=(w(1:1:501))';%取其前一半，并化為列向量mag=abs(H);%求其幅特性db=20*log10((mag+eps)/max(mag));%化為分貝值pha=angle(H);%求其相特性grd=grpdelay(bd,ad,w);%求其群延遲特性subplot(2,2,1);plot(w/pi,mag);title('幅度響應')ylabel('|H|');axis([0,1,0,1.1]);set(gca,'XTickMode','manual','XTick',[0,0.2,0.4,1]);gridon%顯示刻度線set(gca,'YTickMode','manual','YTick',[0,Attn,Ripple,1]);gridonsubplot(2,2,3);plot(w/pi,db);title('幅度(單位：dB)')xlabel('頻率(單位：pi)');ylabel('分貝');axis([0,1,-100,50])set(gca,'XTickMode','manual','XTick',[0,0.2,0.4,1]);gridonset(gca,'YTickMode','manual','YTick',[-100,-50,0,50]);gridon%set(gca,'YTickLabelMode','manual','YTickLabels',['-100';'-50';'0';'50']);subplot(2,2,2);plot(w/pi,pha/pi);title('相位響應')xlabel('');ylabel('單位：pi');axis([0,1,-1,1])set(gca,'XTickMode','manual','XTick',[0,0.2,0.4,1]);gridonset(gca,'YTickMode','manual','YTick',[-1,0,1]);gridonsubplot(2,2,4);plot(w/pi,grd);title('群延遲')xlabel('頻率(單位：pi)');ylabel('樣本');axis([0,1,0,10])set(gca,'XTickMode','manual','XTick',[0,0.2,0.4,1]);gridonset(gca,'YTickMode','manual','YTick',[0:1:10]);gridonset(gcf,'color','w')%置圖形背景色為白色調(diào)試程序；答：運行結(jié)果如下圖3.1至3.2:圖3.1圖3.2要求設計數(shù)字巴特沃思低通濾波器，給定抽樣頻率，要求在頻率小于的通帶內(nèi)，幅度特性下降小于；在頻率大于的阻帶內(nèi)，衰減大于答：源程序如下：clearall;fp=1000;fs=1500;Fs=10000;rp=1;rs=15;%wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;Fs=Fs/Fs;%Firstlytofinishfrequencyprewarping;wap=tan(wp/2);was=tan(ws/2);[n,wn]=buttord(wap,was,rp,rs,'s')%Note:'s'![z,p,k]=buttap(n);[bp,ap]=zp2tf(z,p,k)[bs,as]=lp2lp(bp,ap,wap)%Note:s=(2/Ts)(z-1)(z+1);Ts=1,thatis2fs=1,fs=0.5;[bz,az]=bilinear(bs,as,Fs/2)[h,w]=freqz(bz,az,1256,Fs*10000);plot(w,abs(h));gridon;輸出設計結(jié)果，以及對應的幅頻特性圖并記錄；答：運行結(jié)果如下圖3.3至3.5:圖3.3圖3.4圖3.5四、思考如何編寫以巴特沃思模擬低通濾波器為樣本濾波器，設計基于雙線性變換法的低通數(shù)字帶通濾波器？如何直接利用matlab工具箱設計數(shù)字濾波器？實驗四FIR數(shù)字濾波器的設計實驗目的掌握FIR數(shù)字濾波器設計的方法；熟悉軟件編程實現(xiàn)的方法；認識到FIR數(shù)字濾波器的線性相位的概念；了解IIR和FIR數(shù)字濾波器的各自特點；實驗要求了解和掌握用窗函數(shù)法設計FIR數(shù)字濾波器的過程；整理設計算法和設計步驟；編寫和調(diào)試程序；實驗步驟利用matlab工具箱產(chǎn)生各種窗函數(shù)；（1）矩形窗（RectangleWindow）調(diào)用格式：w=boxcar(n)，根據(jù)長度n產(chǎn)生一個矩形窗w如圖4.1。n=51;wvtool(boxcar(n))圖4.1（2）三角窗（TriangularWindow）調(diào)用格式：w=triang(n)，根據(jù)長度n產(chǎn)生一個三角窗w。n=51;wvtool(triang(n))圖形如下圖4.2：圖4.2（3）漢寧窗（HanningWindow）調(diào)用格式：w=hanning(n)，根據(jù)長度n產(chǎn)生一個漢寧窗w如圖4.3。n=51;wvtool(hanning(n))圖4.3（4）海明窗（HammingWindow）調(diào)用格式：w=hamming(n)，根據(jù)長度n產(chǎn)生一個海明窗w如圖4.4。n=51;wvtool(hamming(n))圖4.4（5）布拉克曼窗（BlackmanWindow）調(diào)用格式：w=blackman(n)，根據(jù)長度n產(chǎn)生一個布拉克曼窗w如圖4.5。程序如下：n=51;wvtool(Blackman(n))圖4.5（6）愷撒窗（KaiserWindow）調(diào)用格式：w=kaiser(n,beta)，根據(jù)長度n和影響窗函數(shù)旁瓣的β參數(shù)產(chǎn)生一個愷撒窗w如圖4.6。程序如下：n=51;beta=8;wvtool(Kaiser(n,beta))圖4.62.編寫用窗函數(shù)法設計FIR低通濾波器程序，調(diào)試；答：源程序如下：wp=0.2*pi;ws=0.4*pi;deltaw=ws-wp;%確定過度帶寬N0=ceil(6.6*pi/deltaw)+1;%確定濾波器階數(shù)N0，ceil是向上取整函數(shù)N=N0+mod(N0+1,2);wdham=(hamming(N))';wc=(ws+wp)/2;%理想低通的截止頻率tao=(N-1)/2;n=[0:N-1];m=n-tao+eps;hd=sin(wc*m)./(pi*m);%求理想脈沖響應h=hd.*wdham;[H,w]=freqz(h,[1],1000,'whole');H=(H(1:1:501))';w=(w(1:1:501))';mag=abs(H);db=20*log10((mag+eps)./max(mag));pha=angle(H);grd=grpdelay(h,[1],w);dw=2*pi/1000;Rp=-(min(db(1:wp/dw+1)));%實際通帶波動As=-round(max(db(ws/dw+1:501)));%最小阻帶衰減stem(n,hd,'.');title('理想沖激響應')axis([0N-1-0.20.3]),ylabel('hd(n)');text(N+1,-0.1,'n')%畫圖subplot(2,2,1);stem(n,hd,'.');title('理想沖激響應')axis([0N-1-0.20.3]),ylabel('hd(n)');text(N+1,-0.1,'n')subplot(2,2,2);stem(n,wdham,'.');title('海明窗')axis([0N-101.1]);ylabel('w(n)');text(N+1,0,'n')subplot(2,2,3);stem(n,h,'.');title('實際沖激響應')axis([0N-1-0.2,0.3]);xlabel('n');ylabel('h(n)')subplot(2,2,4);plot(w/pi,db);title('幅度響應(單位：dB)');gridaxis([01-10010]);xlabel('頻率(單位：pi)');ylabel('分貝數(shù)')n=0;N-1;得到的結(jié)果如圖4.7：圖4.73.設計指標為：抽樣頻率為，通帶截止頻率為，阻帶起始頻率為，阻帶衰減不小于－50dB；答：源程序如下：clearall;Wp=0.2*pi;Ws=0.4*pi;tr_width=Ws-Wp;%過渡帶寬度N=ceil(6.6*pi/tr_width)+1%濾波器長度n=0:1:N-1;Wc=(Ws+Wp)/2;%理想低通濾波器的截止頻率hd=ideal_lp1(Wc,N);%理想低通濾波器的單位沖激響應w_ham=(hamming(N))';%海明窗h=hd.*w_ham;%截取得到實際的單位脈沖響應[db,mag,pha,w]=freqz_m2(h,[1]);%計算實際濾波器的幅度響應delta_w=2*pi/1000;Ap=-(min(db(1:1:Wp/delta_w+1)))%實際通帶紋波As=-round(max(db(Ws/delta_w+1:1:501)))%實際阻帶紋波subplot(221)stem(n,hd)title('理想單位脈沖響應hd(n)')subplot(222)stem(n,w_ham)title('海明窗w(n)')subplot(223)stem(n,h)title('實際單位脈沖響應hd(n)')subplot(224)plot(w/pi,db)title('幅度響應(dB)')axis([0,1,-100,10])其中freqz_m2函數(shù)定義如下：function[db,mag,pha,w]=freqz_m2(b,a)%濾波器的幅值響應(相對、絕對)、相位響應%db:相對幅值響應%mag：絕對幅值響應%pha:相位響應%w采樣頻率；%b系統(tǒng)函數(shù)H(z)的分子項(對FIR，b=h)%a系統(tǒng)函數(shù)H(z)的分母項(對FIR，a=1)[H,w]=freqz(b,a,1000,'whole');H=(H(1:1:501))';w=(w(1:1:501))';mag=abs(H);%絕對幅值響應db=20*log10((mag+eps)/max(mag));%相對幅值響應pha=angle(H);%相位響應另外還有函數(shù)ideal_lp1：functionhd=ideal_lp1(Wc,N)%computetheideallowpassfiterunitpulserespondencehd(n)%wc:cutofffrequency%N:windowlength%hd:unitpulserespondencealpha=(N-1)/2;n=0:1:N-1;m=n-alpha+eps;hd=sin(Wc*m)./(pi*m);輸出設計結(jié)果，繪制對應的幅頻特性圖，并記錄；答：得到結(jié)果如下圖4.8：圖4.8思考如何用窗函數(shù)法設計其他類型的線性相位的FIR濾波器，如高通、帶通和帶阻濾波器？答：這個問題將在實驗五的對比過程中有所體現(xiàn)。實驗五數(shù)字信號處理綜合設計一、實驗目的1．學會MATLAB的使用，掌握MATLAB的程序設計方法；2．掌握在Windows環(huán)境下語音信號采集的方法；3．掌握數(shù)字信號處理的基本概念、基本理論和基本方法；4．掌握MATLAB設計FIR和IIR數(shù)字濾波器的方法；5．學會用MATLAB對信號進行分析和處理。二、實驗要求和步驟1．語音信號的采集要求利用windows下的錄音機或其他軟件，錄制一段自己的話音，時間控制在1秒左右。然后在MATLAB軟件平臺下，利用函數(shù)wavread對語音信號進行采樣，記住采樣頻率和采樣點數(shù)。通過wavread函數(shù)的使用，要求理解采樣頻率、采樣位數(shù)等概念。wavread函數(shù)調(diào)用格式：y=wavread(file)，讀取file所規(guī)定的wav文件，返回采樣值放在向量y中。[y,fs,nbits]=wavread(file)，采樣值放在向量y中，fs表示采樣頻率（Hz），nbits表示采樣位數(shù)。y=wavread(file,N)，讀取前N點的采樣值放在向量y中。y=wavread(file,[N1,N2])，讀取從N1點到N2點的采樣值放在向量y中。答：信號采集源程序：[y,fs,bits]=wavread('方品語音.wav');sound(y,fs,bits);Y=fft(y,4096);figure(1)plot(y)title('原時域波形');ylabel('幅度');xlabel('n');figure(2)plot(abs(Y))axis([0,4096,0,3]);title('原頻譜特性');ylabel('幅度');xlabel('頻率/hz');采集到的時域信號和頻譜特性如下圖5.1至5.2：圖5.1時域波形圖5.2頻譜特性2．語音信號的頻譜分析要求首先畫出語音信號的時域波形；然后對語音信號進行頻譜分析，在MATLAB中，可以利用函數(shù)fft對信號進行快速付立葉變換，得到信號的頻譜特性；從而加深對頻譜特性的理解。答：如圖5.2.3．設計數(shù)字濾波器和畫出頻率響應根據(jù)語音信號的特點給出有關濾波器的性能指標：1）低通濾波器性能指標，fp=1000Hz，fc=1200Hz，As=100dB，Ap=1dB；2）高通濾波器性能指標，fc=4800Hz，fp=5000HzAs=100dB，Ap=1dB；3）帶通濾波器性能指標，fp1=1200Hz，fp2=3000Hz，fc1=1000Hz，fc2=3200Hz，As=100dB，Ap=1dB。學生可以用窗函數(shù)法設計上面要求的三種濾波器，在MATLAB中，可以利用函數(shù)fir1設計FIR濾波器；或者也可以用雙線性變換法設計上面要求的三種濾波器，在MATLAB中，可以利用函數(shù)butte、cheby1和ellip設計IIR濾波器；最后，利用MATLAB中的函數(shù)freqz畫出各濾波器的頻率響應。答：1、利用窗函數(shù)法設計這三種濾波器：由于所給濾波器的模擬指標要求阻帶的最小衰減As=100dB,所以只能選擇kaiser窗，且kaiser窗的參數(shù)a=10.056，過渡帶寬w為10.8*pi，選取階數(shù)N=49。窗函數(shù)法低通濾波器設計程序及頻率響應如下:①源程序如下：wn=kaiser(49);fc=1200;fs=8000;wc=2*fc/fs;b=fir1(48,wc,wn);freqz(b,1)freqz(b);title('窗函數(shù)低通濾波器頻率響應')②獲得頻率響應如下圖5.3：圖5.3窗函數(shù)設計低通濾波器窗函數(shù)法帶通濾波器設計程序及頻率響應如下:①源程序如下：wn=kaiser(49);fc1=1000;fc2=3200;fs=8000;wc1=2*fc1/fs;wc2=2*fc2/fs;b=fir1(48,[wc1wc2],wn);freqz(b,1)title('窗函數(shù)帶通濾波器響應')②獲得頻率響應如下圖5.4：圖5.4窗函數(shù)設計帶通濾波器窗函數(shù)法高通濾波器設計程序及頻率響應如下:①源程序如下：wn=kaiser(49);fc=4800;wc=2*fc2/fs;fs=22050;b=fir1(48,wc,'high',wn);freqz(b,1);figure(1);freqz(b);title('窗函數(shù)高通濾波器響應')②獲得頻率響應如下圖5.5：圖5.5窗函數(shù)設計高通濾波器利用雙線性法設計這三種濾波器：雙線性變換法設計巴特沃斯低通濾波器，程序和頻率響應如下：①源程序如下：FS=8000;F=1000;Fh=1200;wp=(F/FS)*2*pi;ws=(Fh/FS)*2*pi;OmegaP=2*FS*tan(wp/2);OmegaS=2*FS*tan(ws/2);[n,Wn]=buttord(OmegaP,OmegaS,1,100,'s');[b,a]=butter(n,Wn,'s');[bz,az]=bilinear(b,a,FS);freqz(bz,az,4096,FS,'whole');title('雙線性變換法低通濾波器')②獲得頻率響應如下圖5.6：圖5.6雙線性變換法的低通濾波器頻譜響應（2）雙線性變換法設計巴特沃斯高通濾波器，程序和頻率響應如下：①源程序如下：FS=8000;Fl=2800;Fh=3000;wp=(Fh/FS)*2*pi;ws=(Fl/FS)*2*pi;OmegaP=2*FS*tan(wp/2);OmegaS=2*FS*tan(ws/2);[n,Wn]=buttord(OmegaP,OmegaS,1,100,'s');[b,a]=butter(n,Wn,'s');[bz,az]=bilinear(b,a,FS);freqz(bz,az,4096,FS,'whole');title('雙線性變換法高通濾波器')②獲得頻率響應如下圖5.7：圖5.7雙線性變換法高通濾波器雙線性變換法設計切比雪夫帶通濾波器，程序和頻率響應如下：①源程序如下：wp1=1200/8000*2*pi;wp2=3000/8000*2*pi;ws1=1000/8000*2*pi;ws2=3200/8000*2*pi;Rp=1;Rs=100;Wp1=tan(wp1/2);Wp2=tan(wp2/2);Ws1=tan(ws1/2);Ws2=tan(ws2/2);BW=Wp2-Wp1;W0=Wp1*Wp2;W00=sqrt(W0);WP=1;WS=WP*(W0^2-Ws1^2)/(Ws1*BW);[N,Wn]=cheb1ord(WP,WS,Rp,Rs,'s');[B,A]=cheby1(N,Rp,Wn,'s');[BT,AT]=lp2bp(B,A,W00,BW);[num,den]=bilinear(BT,AT,0.5);[h,omega]=freqz(num,den,32);subplot(2,1,1);stem(omega/pi,abs(h));xlabel('\omega/\pi')ylabel('|H(z)|');subplot(2,1,2);stem(omega/pi,20*log10(abs(h)));xlabel('\omega/\pi')ylabel('增益db');②獲得頻率響應如下圖5.8：圖5.8雙線性變換法帶通濾波器4．用濾波器對信號進行濾波比較兩種濾波器的性能，然后用性能好的各濾波器分別對采集的信號進行濾波，在MATLAB中，F(xiàn)IR濾波器利用函數(shù)fftfilt對信號進行濾波，IIR濾波器利用函數(shù)filter對信號進行濾波。答：1、濾波程序（1）巴特沃斯低通濾波器濾波程序：①源程序如下：[y,fs,bits]=wavread('方品語音.wav');sound(y,fs,bits);Y=fft(y,4096);FS=8000;F=1000;Fh=1200;wp=(F/FS)*2*pi;ws=(Fh/FS)*2*pi;OmegaP=2*FS*tan(wp/2);OmegaS=2*FS*tan(ws/2);[n,Wn]=buttord(OmegaP,OmegaS,1,100,'s');[b,a]=butter(n,Wn,'s');[bt,at]=bilinear(b,a,FS);%用雙線性變換法將模擬轉(zhuǎn)換成數(shù)字帶通濾波器freqz(bt,at,4096,FS,'whole');x=filter(b,a,y);X=fft(x,4096);figure(2);subplot(2,1,1);plot(abs(Y));axis([0,4096,0,3]);title('雙線性巴特沃斯低通濾波器前頻譜');subplot(2,1,2);plot(abs(X));axis([0,4096,0,3]);title('雙線性巴特沃斯低通濾波器后頻譜');sound(x,fs,bits)②獲得頻率響應如下圖5.9：如果程序正常應該如圖5.10：圖5.10理想要出的結(jié)果巴特沃斯高通濾波器濾波程序①源程序如下：[y,fs,bits]=wavread('方品語音.wav');sound(y,fs,bits);Y=fft(y,4096);FS=8000;Fl=2800;Fh=3200;wp=(Fh/FS)*2*pi;ws=(Fl/FS)*2*pi;OmegaP=2*FS*tan(wp/2);OmegaS=2*FS*tan(ws/2);[n,Wn]=buttord(OmegaP,OmegaS,1,100,'s');[b,a]=butter(n,Wn,'s');[bz,az]=bilinear(b,a,FS);freqz(bz,az,4096,FS,'whole');x=filter(b,a,y);X=fft(x,4096);figure(2);subplot(2,1,1);plot(abs(Y));axis([0,4096,0,3]);title('雙線性巴特沃斯高通濾波器前頻譜');subplot(2,1,2);plot(abs(X));axis([0,4096,0,3]);title('雙線性巴特沃斯高通濾波器后頻譜');sound(x,fs,bits)②獲得頻率響應如下圖5.11：圖5.11如果程序正常應該如圖5.12：圖5.12程序正常出的波形分析：以上兩個程序出現(xiàn)警告，可是我通過修改參數(shù)還是不能避免警告，不知道程序問題出在哪里？切比雪夫帶通濾波器濾波程序：①源程序如下：[y,fs,bits]=wavread('方品語音.wav');sound(y,fs,bits);Y=fft(y,4096);FS=8000;Fl=2800;Fh=3200;wp=(Fh/FS)*2*pi;ws=(Fl/FS)*2*pi;OmegaP=2*FS*tan(wp/2);OmegaS=2*FS*tan(ws/2);[n,Wn]=buttord(OmegaP,OmegaS,1,100,'s');[b,a]=butter(n,Wn,'s');[bz,az]=bilinear(b,a,FS);freqz(bz,az,4096,FS,'whole');x=filter(b,a,y);X=fft(x,4096);figure(2);subplot(2,1,1);plot(abs(Y));axis([0,4096,0,3]);title('雙線性巴特沃斯高通濾波器前頻譜');subplot(2,1,2);plot(abs(X));axis([0,4096,0,3]);title('雙線性巴特沃斯高通濾波器后頻譜');sound(x,fs,bits)②獲得頻率響應如下圖5.13：圖5.13（4）窗函數(shù)設計低通濾波器濾波程序:①源程序如下：[y,fs,bits]=wavread('方品語音.wav');sound(y,fs,bits);Y=fft(y,4096);wn=kaiser(49);fc=1200;fs=8000;wc=2*fc/fs;b=fir1(48,wc,wn);freqz(b,1)x=fftfilt(b,y);X=fft(x,4096);subplot(2,1,1);plot(abs(Y));axis([0,4096,0,3]);title('窗函數(shù)低通濾波器前頻譜');subplot(2,1,2);plot(abs(X));axis([0,4096,0,3]);title('窗函數(shù)低通濾波器后頻譜');sound(x,fs,bits)figure(2);freqz(b);title('窗函數(shù)低通濾波器頻率響應')figure(3);subplot(2,1,1);plot(abs(y));axis([0,4096,0,1.2]);title('窗函數(shù)低通濾波器前波形');subplot(2,1,2);plot(abs(x));axis([0,4096,0,1.2]);title('窗函數(shù)低通濾波器后波形');sound(x,fs,bits)②獲得頻率響應如下圖5.14至5.16：圖5.14圖5.15圖5.16（5）窗函數(shù)高通濾波器濾波程序：①源程序如下：wn=kaiser(49);fc=4800;wc=2*fc2/fs;fs=22050;b=fir1(48,wc,'high',wn);freqz(b,1);x=fftfilt(b,y);X=fft(x,4096);subplot(2,1,1);plot(abs(y));axis([0,4096,0,1]);title('窗函數(shù)高通濾波器前波形');subplot(2,1,2);plot(abs(x));axis([0,4096,0,0.3]);title('窗函數(shù)gao通濾波器后波形');sound(x,fs,bits)figure(2);subplot(2,1,1);plot(abs(Y));axis([0,4096,0,3]);title('窗函數(shù)高通濾波器前頻譜');subplot(2,1,2);plot(abs(X));axis([0,4096,0,2.5]);title('窗函數(shù)高通濾波器后頻譜');sound(x,fs,bits)figure(3)freqz(b);title('窗函數(shù)高通濾波器頻率響應')②獲得頻率響應如下圖5.17至5.19：圖5.17圖5.18圖5.19（6）窗函數(shù)設計帶通濾波器程序①源程序如下：[y,fs,bits]=wavread('方品語音.wav');sound(y,fs,bits);Y=fft(y,4096);wn=kaiser(49);fc1=1000;fc2=3200;fs=8000;wc1=2*fc1/fs;wc2=2*fc2/fs;b=fir1(48,[wc1wc2],wn);freqz(b,1)x=fftfilt(b,y);X=fft(x,4096);subplot(2,1,1);plot(abs(y));axis([0,4096,0,1]);title('窗函數(shù)帶通濾波器前波形');subplot(2,1,2);plot(abs(x));axis([0,4096,0,1]);title('窗函數(shù)帶通濾波器后波形');sound(x,fs,bits)figure(2);subplot(2,1,1);plot(abs(Y));axis([0,4096,0,3]);title('窗函數(shù)帶通濾波器前頻譜');subplot(2,1,2);plot(abs(X));axis([0,4096,0,3]);title('窗函數(shù)帶通濾波器后頻譜');sound(x,fs,bits)figure(3);freqz(b);title('窗函數(shù)帶通濾波器頻率響應')②獲得頻率響應如下圖5.20至5.22：圖5.20圖5.21圖5.22分析：對比（1）和（4）即雙線性法巴特沃斯低通濾波器和窗函數(shù)法低通濾波器對比可發(fā)現(xiàn)：窗函數(shù)濾波后頻譜阻帶的波紋較大一些，而巴特沃斯阻帶較平滑；對比（2）和（5）即巴特沃斯高通和窗函數(shù)法高通濾波器對比可發(fā)現(xiàn)：窗函數(shù)在阻帶的波紋較大，即濾波效果不如巴特沃斯濾波器；對比（3）和（6）的圖可知切比雪夫濾波頻譜比較尖銳，只能通過很小頻段的波形，而窗函數(shù)能通過適當范圍的頻率波形，且濾波效果較為理想。5．比較濾波前后語音信號的波形及頻譜要求在一個窗口同時畫出濾波前后的波形及頻譜。答：對比結(jié)果如上文。6．回放語音信號在MATLAB中，函數(shù)sound可以對聲音進行回放。其調(diào)用格式：sound(x,fs,bits)；可以感覺濾波前后的聲音有變化。三、實驗報告要求1.簡述實驗原理及目的。2.按照實驗步驟及要求，比較各種情況下的濾波性能。3.總結(jié)實驗所得主要結(jié)論。4.簡要回答思考題。四、思考1．雙線性變換法中Ω和ω之間的關系是非線性的，在實驗中你注意到這種非線性關系了嗎？從那幾種數(shù)字濾波器的幅頻特性曲線中可以觀察到這種非線性關系？答：注意到了。雙線性變換消除了混疊，但同時帶來了嚴重的非線性失真。從雙線性變換法設計各種濾波器的幅頻特性中均可看住這種非線性關系。2．能否利用公式完成脈沖響應不變法的數(shù)字濾波器設計？為什么？答：能。脈沖響應不變法的設計原理是利用數(shù)字濾波器的單位抽樣響應序列H(z)來逼近模擬濾波器的沖激響應g(t)。按照脈沖響應不變法的原理，通過模擬濾波器的系統(tǒng)傳遞函數(shù)G(s)，可以直接求得數(shù)字濾波器的系統(tǒng)函數(shù)H(z)，其轉(zhuǎn)換步驟如下：1)利用ω=ΩT(可由關系式推導出)，將，轉(zhuǎn)換成，而不變；2)求解低通模擬濾波器的傳遞函數(shù)G(s)；3)將模擬濾波器的傳遞函數(shù)G(s)轉(zhuǎn)換為數(shù)字濾波器的傳遞函數(shù)H(z)。實驗６Simulink部分一、實驗目的1．學會MATLAB中Simulink的使用；2．掌握在Simulink的信號產(chǎn)生方法；3．掌握Simulink中的離散傅立葉變換的計