數(shù)字信號(hào)實(shí)驗(yàn)報(bào)告1

數(shù)字信號(hào)處理導(dǎo)論實(shí)驗(yàn)報(bào)告姓名：金濤 學(xué)號(hào)：201005090209實(shí)驗(yàn)一信號(hào)、系統(tǒng)及系統(tǒng)響應(yīng)當(dāng)系統(tǒng)的輸入輸出差分方程為：Y(n)-0.8y(n-1)-0.5y(n-2)=0.7x(n)+0.3x（n-1）通過MATLAB編程實(shí)現(xiàn)并考慮如下問題：（1）當(dāng)系統(tǒng)的輸入為單位沖激函數(shù)時(shí)，分別利用filter函數(shù)和impz函數(shù)得到的系統(tǒng)單位沖激響應(yīng)曲線。（2）當(dāng)系統(tǒng)的輸入為單位階躍函數(shù)時(shí)，分別利用filter函數(shù)和impz函數(shù)得到的系統(tǒng)單位階躍響應(yīng)曲線。（3）對(duì)于不同的輸入，系統(tǒng)的輸出有什么變化，試討論之。實(shí)驗(yàn)原理對(duì)一個(gè)給定的LSI系統(tǒng)，其轉(zhuǎn)移函數(shù)H(z)的定義和表示形式為：.習(xí)慣上，令H(z)=B(z)/A(z).在MATLAB中，因?yàn)閿?shù)組的下標(biāo)不能為零（當(dāng)然也不能為負(fù)值），因此，可重新表示為H(z)=在有關(guān)MATLAB的系統(tǒng)分析的文件中。分子和分母的系數(shù)被定義為向量，即b=[b(1),b(2),b(3),...,b(+1)]a=[a(1),a(2),...,a()]并要求a（1）=1，如果a(1)≠1,則程序自動(dòng)將其規(guī)劃為1實(shí)驗(yàn)內(nèi)容源程序x=ones(100);t=1:100;%產(chǎn)生單位階躍序列b=[.7,.3];%b向量a=[1,.8,.5];%a向量y=filter(b,a,x);%實(shí)現(xiàn)實(shí)驗(yàn)1，圖（1）plot(t,x,'g.',t,y,'k-');[h,t]=impz(b,a,40);%求出單位抽樣響應(yīng)。圖（2）stem(t,h,'.');gridon;t=0:20;x=[1,zeros(1,20)];%產(chǎn)生單位采樣序列b=[.7,.3];%形成b向量a=[1,.8,.5];%形成a向量y=filter(b,a,x);%filter函數(shù)圖（3）plot(t,x,'g.',t,y,'k-');t=0:20;x=[1,zeros(1,20)];b=[.7,.3];a=[1,.8,.5];[h,t]=impz(b,a,40);%impz函數(shù)圖（4）stem(t,h,'.');gridon;實(shí)驗(yàn)結(jié)果（1）（2）（3）（4）分析結(jié)果輸入為離散是，輸出為連續(xù)。相反，輸入為連續(xù)，輸出為離散。對(duì)于單位階躍和單位抽樣輸入來說，輸出沒有變化實(shí)驗(yàn)二使用FFT作頻譜分析（1）使用FFT對(duì)MATLAB中的數(shù)據(jù)noissin.dat進(jìn)行譜分析。（2）使用FFT對(duì)語音數(shù)據(jù)noisyspeech.wav進(jìn)行譜分析。實(shí)驗(yàn)原理（1）離散傅里葉變換(DFT)公式為：X(k)=∑x(n)W^nk;x(n)=∑X(k)W^-nk;其中w=e^(2∏nk/N),N為離散序列的長(zhǎng)度。（2）快速傅里葉變換(FFT)是利用w因子的取值特點(diǎn)，減少DFT的復(fù)數(shù)乘法的次數(shù)。其中一種是時(shí)間抽取基2算法，它將時(shí)間按奇偶逐級(jí)分開，直到兩點(diǎn)的DFT。MATLAB提供了fft函數(shù)可用于計(jì)算FFT，器調(diào)用形式為;X=fft（x）或X=fft（x，N），N為2的整數(shù)次冪，若x的長(zhǎng)度小于N，則補(bǔ)零，若超過則舍去N后的數(shù)據(jù)。（3）函數(shù)形式[y,fs,bits]=wavread('Blip',[N1N2]);用于讀取語音，采樣值放在向量y中，fs表示采樣頻率(Hz)，bits表示采樣位數(shù)。[N1N2]表示讀取從N1點(diǎn)到N2點(diǎn)的值（若只有一個(gè)N的點(diǎn)則表示讀取前N點(diǎn)的采樣值）。函數(shù)形式s=noissin（n1，n2）用于讀取MALAB的噪聲信號(hào)。實(shí)驗(yàn)內(nèi)容xx=wavread('noisyspeech.wav');fs=100;N=128;x=xx(1:N);n=1:N;X=abs(fft(x,N));subplot(221);plot(n,x);xlabel('n');ylabel('x(n)');gridon;subplot(222);plot(n,X);xlabel('k');ylabel('|X(k)|');gridon;loadnoissin;s=noissin(1:20);S=fft(s);subplot(223);plot(abs(s));xlabel('n');ylabel('|s(n)|');gridon;subplot(224);plot(abs(S));xlabel('k');ylabel('|S(k)|');gridon;實(shí)驗(yàn)結(jié)果及分析實(shí)驗(yàn)三使用雙線性Z變換設(shè)計(jì)IIR濾波器使用雙線性Z變換法設(shè)計(jì)一個(gè)低通數(shù)字IIR濾波器，給定的數(shù)字濾波器的技術(shù)指標(biāo)為fp=100Hz,fs=300Hz,ap=3dB,as=20dB,抽樣頻率Fs=1000Hz實(shí)驗(yàn)原理1）設(shè)計(jì)濾波器就是要設(shè)計(jì)一個(gè)系統(tǒng)是其能讓一定頻率的波段通過或?yàn)V去，對(duì)IIR濾波器，器轉(zhuǎn)移函數(shù)是：H(Z)=(∑bz^(-r))/(1+∑az^(-k))。（2）設(shè)計(jì)的一般原則：若使濾波器拒絕某個(gè)頻率，應(yīng)在單位園上相應(yīng)的頻率處設(shè)置一個(gè)零點(diǎn)，反之則設(shè)置一個(gè)極點(diǎn)。(3)低通數(shù)字IIR濾波器設(shè)計(jì)步驟：給出數(shù)字低通濾波器的技術(shù)要求→映射為模擬低通的技術(shù)要求→歸一化為模擬低通濾波器的技術(shù)要求→設(shè)計(jì)出G(P)→G(S)→映射到數(shù)字濾波器的轉(zhuǎn)移函數(shù)G(Z)。(3)雙線性Z變換，即S平面到Z平面的映射關(guān)系：S=2（Z-1）/Ts（Z+1）。實(shí)驗(yàn)內(nèi)容fp=100;fs=300;Fs=1000;%給定要設(shè)計(jì)的低通濾波器的設(shè)計(jì)頻率；rp=3;rs=20;%給定需要的衰減，單位為db；wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;Fs=Fs/Fs;%令Fs=1wap=tan(wp/2);was=tan(ws/2);%實(shí)現(xiàn)由數(shù)字濾波器到模擬濾波器的轉(zhuǎn)換；[n,wn]=buttord(wap,was,rp,rs,'s');%求模擬濾波器的階次；[z,p,k]=buttap(n);%設(shè)計(jì)模擬低通巴特沃斯濾波器。得到極零點(diǎn)；[bp,ap]=zp2tf(z,p,k);%由極零點(diǎn)得到濾波器分子分母多項(xiàng)式的系數(shù)；[bs,as]=lp2lp(bp,ap,wap);%歸一化低通到實(shí)際低通的轉(zhuǎn)換；[bz,az]=bilinear(bs,as,Fs/2);%實(shí)現(xiàn)s到z的轉(zhuǎn)換；[h,w]=freqz(bz,az,256,Fs*1000);%求出濾波器的頻率響應(yīng)；plot(w,abs(h));gridon;實(shí)驗(yàn)結(jié)果結(jié)果分析由圖可知，數(shù)字濾波器的幅頻曲線完全滿足技術(shù)要求，而且在f>fp后，數(shù)字濾波器的幅頻響應(yīng)的衰減比較快，這正是我們希望的實(shí)驗(yàn)四窗函數(shù)法設(shè)計(jì)FIR數(shù)字濾波器根據(jù)下列指標(biāo)，設(shè)計(jì)一個(gè)FIR數(shù)字低通濾波器：ω=2.0∏，ω=4.0∏Ap=0.25dB，As=50dB。（1）分別考查矩形窗、三角窗、Hanning窗、Hamming窗并分析這些不同的窗函數(shù)對(duì)濾波器的幅度響應(yīng)有什么影響。（2）選擇一個(gè)合適的窗函數(shù)，確定單位沖激響應(yīng)，繪出所設(shè)計(jì)的濾波器的幅度響應(yīng)。（一）實(shí)驗(yàn)?zāi)康模毫私獬Ｓ玫膸追N窗函數(shù)，能正確選擇適當(dāng)?shù)拇昂瘮?shù)進(jìn)行濾波器設(shè)計(jì)；掌握窗函數(shù)法設(shè)計(jì)數(shù)字低通濾波器。（二）實(shí)驗(yàn)原理常用的窗函數(shù)：1、矩形窗函數(shù)為boxcar和rectwin，調(diào)用格式：w=boxcar（N）w=rectwin(N)其中N是窗函數(shù)的長(zhǎng)度，返回值w是一個(gè)N階的向量。2、三角窗函數(shù)為triang，調(diào)用格式：w=triang(N)3、漢寧窗函數(shù)為hann，調(diào)用格式：w=hann（N）4、海明窗函數(shù)為hamming，調(diào)用格式：w=hamming（N）各個(gè)窗函數(shù)的性能比較窗函數(shù)第一瓣相對(duì)于主瓣衰減（dB）主瓣寬阻帶最小衰減矩形窗-13

4∏/N21三角窗-25

8∏/N25漢寧窗-31

8∏/N44海明窗-41

8∏/N53（三）實(shí)驗(yàn)內(nèi)容1、題一：生成四種窗函數(shù)。矩形窗、三角窗、漢寧窗、海明窗，并觀察其頻率響應(yīng)。源程序：clearall;n=20;window=rectwin(n);[h,w]=freqz(window,1);subplot(4,2,1);stem(window);title('矩形窗');subplot(4,2,2);plot(w/pi,20*log(abs(h))/abs(h(1)));title('矩形窗的頻率響應(yīng)');holdon;window=triang(n);[h,w]=freqz(window,1);subplot(4,2,3);stem(window);title('三角窗');subplot(4,2,4);plot(w/pi,20*log(abs(h))/abs(h(1)));title('三角窗的頻率響應(yīng)');holdon;window=hann(n);[h,w]=freqz(window,1);subplot(4,2,5);stem(window);title('漢寧窗');subplot(4,2,6);plot(w/pi,20*log(abs(h))/abs(h(1)));title('漢寧窗的頻率響應(yīng)');holdon;window=hamming(n);[h,w]=freqz(window,1);subplot(4,2,7);stem(window);title('海明窗');subplot(4,2,8);plot(w/pi,20*log(abs(h))/abs(h(1)));title('海明窗的頻率響應(yīng)');holdon;實(shí)驗(yàn)結(jié)果：分析實(shí)驗(yàn)結(jié)果：通過對(duì)窗函數(shù)的調(diào)用，使4個(gè)窗函數(shù)的時(shí)域及頻域圖放置在一起，比較直觀地看出各個(gè)窗函數(shù)的特性。2、題二：根據(jù)下列技術(shù)指標(biāo)，設(shè)計(jì)一個(gè)FIR數(shù)字低通濾波器：wp=0.2π，ws=0.4π，ap=0.25dB，as=50dB，選擇一個(gè)適當(dāng)?shù)拇昂瘮?shù)，確定單位沖激響應(yīng)，繪出所設(shè)計(jì)的濾波器的幅度響應(yīng)。提示：根據(jù)窗函數(shù)最小阻帶衰減的特性表，可采用海明窗可提供大于50dB的衰減，其過渡帶為6.6π/N，因此具有較小的階次。源程序：clearall;wp=0.2*pi;%通帶截止頻率ws=0.4*pi;%阻帶截止頻率tr_wdith=ws-wp;%過渡帶寬度N=ceil(6.6*pi/tr_wdith)+1;%得到Nn=0:1:N-1;wc=(ws+wp)/2;%理想低通濾波器的截止頻率hd=ideal_lp1(wc,N);%理想低通濾波器的單位沖激響應(yīng)w_ham=(hamming(N))';%海明窗h=hd.*w_ham;%截取得到實(shí)際的單位脈沖響應(yīng)[db,mag,pha,w]=freqz_m2(h,1);%計(jì)算實(shí)際濾波器的幅度響應(yīng)subplot(221);stem(n,hd,'k');title('理想低通濾波器的單位脈沖響應(yīng)hd（n）');subplot(222);stem(n,w_ham,'k');title('海明窗w（n）');subplot(223);stem(n,h,'k');title('實(shí)際低通濾波器的單位脈沖響應(yīng)h（n）');subplot(224);plot(w/pi,db,'k');title('實(shí)際低通濾波器的幅度響應(yīng)（dB）');axis([0,1,-100,10]);子函數(shù)freqz_m2(b,a)：function[db,mag,pha,w]=freqz_m2(b,a);%濾波器的幅值響應(yīng)（相對(duì)、絕對(duì)）、相位響應(yīng)%db：相對(duì)幅值響應(yīng)；%mag：絕對(duì)幅值響應(yīng)；%pha：相位響應(yīng)；%w：采樣頻率；%b:系統(tǒng)函數(shù)H（Z）的分子項(xiàng)（對(duì)FIR,b=h）;%a:系統(tǒng)函數(shù)H（Z）的分母項(xiàng)（對(duì)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(wc,N)：functionhd=ideal_lp1(wc,N);alpha=(N-1)/2;n=0:1:N-1;m=n-alpha+eps;hd=sin(wc*m)./(pi*m);實(shí)驗(yàn)結(jié)果：分析實(shí)驗(yàn)結(jié)果：用窗函數(shù)的方法設(shè)計(jì)FIR低通數(shù)字濾波器，其性能十分逼近與理想的IIR數(shù)字低通濾波器，并且具有較好的衰減特性。如圖所示。實(shí)驗(yàn)五用C++編程實(shí)現(xiàn)FFT并作頻分析（一）源程序#include<stdio.h>#include<stdlib.h>#include<math.h>constintN=1024;constfloatPI=3.1416;inlinevoidswap(float&a,float&b){floatt;t=a;a=b;b=t;}voidbitrp(floatxreal[],floatximag[],intn){//位反轉(zhuǎn)置換Bit-reversalPermutationinti,j,a,b,p;for(i=1,p=0;i<n;i*=2){p++;}for(i=0;i<n;i++){a=i;b=0;for(j=0;j<p;j++){b=(b<<1)+(a&1);//b=b*2+a%2;a>>=1;//a=a/2;}if(b>i){swap(xreal[i],xreal[b]);swap(ximag[i],ximag[b]);}}}voidFFT(floatxreal[],floatximag[],intn){//快速傅立葉變換，將復(fù)數(shù)x變換后仍保存在x中，xreal,ximag分別是x的實(shí)部和虛部floatwreal[N/2],wimag[N/2],treal,timag,ureal,uimag,arg;intm,k,j,t,index1,index2;bitrp(xreal,ximag,n);//計(jì)算1的前n/2個(gè)n次方根的共軛復(fù)數(shù)W'j=wreal[j]+i*wimag[j],j=0,1,...,n/2-1arg=-2*PI/n;treal=cos(arg);timag=sin(arg);wreal[0]=1.0;wimag[0]=0.0;for(j=1;j<n/2;j++){wreal[j]=wreal[j-1]*treal-wimag[j-1]*timag;wimag[j]=wreal[j-1]*timag+wimag[j-1]*treal;}for(m=2;m<=n;m*=2){for(k=0;k<n;k+=m){for(j=0;j<m/2;j++){index1=k+j;index2=index1+m/2;t=n*j/m;//旋轉(zhuǎn)因子w的實(shí)部在wreal[]中的下標(biāo)為ttreal=wreal[t]*xreal[index2]-wimag[t]*ximag[index2];timag=wreal[t]*ximag[index2]+wimag[t]*xreal[index2];ureal=xreal[index1];uimag=ximag[index1];xreal[index1]=ureal+treal;ximag[index1]=uimag+timag;xreal[index2]=ureal-treal;ximag[index2]=uimag-timag;}}}}voidIFFT(floatxreal[],floatximag[],intn){//快速傅立葉逆變換floatwreal[N/2],wimag[N/2],treal,timag,ureal,uimag,arg;intm,k,j,t,index1,index2;bitrp(xreal,ximag,n);//計(jì)算1的前n/2個(gè)n次方根Wj=wreal[j]+i*wimag[j],j=0,1,...,n/2-1arg=2*PI/n;treal=cos(arg);timag=sin(arg);wreal[0]=1.0;wimag[0]=0.0;for(j=1;j<n/2;j++){wreal[j]=wreal[j-1]*treal-wimag[j-1]*timag;wimag[j]=wreal[j-1]*timag+wimag[j-1]*treal;}for(m=2;m<=n;m*=2){for(k=0;k<n;k+=m){for(j=0;j<m/2;j++){index1=k+j;index2=index1+m/2;t=n*j/m;//旋轉(zhuǎn)因子w的實(shí)部在wreal[]中的下標(biāo)為ttreal=wreal[t]*xreal[index2]-wimag[t]*ximag[index2];timag=wreal[t]*ximag[index2]+wimag[t]*xreal[index2];ureal=xreal[index1];uimag=ximag[index1];xreal[index1]=ureal+treal;ximag[index1]=uimag+timag;xreal[index2]=ureal-treal;ximag[index2]=uimag-timag;}}}for(j=0;j<n;j++){xreal[j]/=n;ximag[j]/=n;}}voidFFT_test(){charinputfile[]="input.txt";//從文件input.txt中讀入原始數(shù)據(jù)charoutputfile[]="output.txt";//將結(jié)果輸出到文件output.txt中floatxreal[N]={},ximag[N]={};intn,i;FILE*input,*output;if(!(input=fopen(inputfile,"r"))){printf("Cannotopenfile.");exit(1);}if(!(output=fopen(outputfile,"w"))){printf("Cannotopenfile.");exit(1);}i=0;while((fscanf(input,"%f%f",xreal+i,ximag+i))!=EOF){i++;}n=i;//要求n為2的整數(shù)冪while(i>1){if(i%2){fprintf(output,"%disnotapowerof2!",n);exit(1);}i/=2;}FFT(xreal,ximag,n);fprintf(output,"FFT:irealimag");for(i=0;i<n;i++){fprintf(output,"%4d%8.4f%8.4f",i,xreal[i],ximag[i]);}fprintf(output,"=================================");IFFT(xreal,ximag,n);fprintf(output,"IFFT:irealimag");for(i=0;i<n;i++){fprintf(output,"%4d%8.4f%8.4f",i,xreal[i],ximag[i]);}if(fclose(input)){printf("Filecloseerror.");exit(1);}if(fclose(output)){printf("Filecloseerror.");exit(1);}}intmain(){FFT_test();return0;}（二）實(shí)