哈工大數(shù)字信號處理報告

word文檔可自由復(fù)制編輯word文檔可自由復(fù)制編輯word文檔可自由復(fù)制編輯數(shù)字信號處理實驗報告課程名稱：數(shù)字信號處理院系：班級：姓名：學(xué)號：指導(dǎo)教師：實驗時間：2013年11月實驗一：用FFT作譜分析一、實驗?zāi)康?．進(jìn)一步加深DFT算法原理和基本性質(zhì)的理解(因為FFT只是DFT的一種快速算法，所以FFT的運算結(jié)果必然滿足DFT的基本性質(zhì))。2．熟悉FFT算法原理和FFT子程序的應(yīng)用。3．學(xué)習(xí)用FFT對連續(xù)信號和時域離散信號進(jìn)行譜分析的方法，了解可能出現(xiàn)的分析誤差及其原因，以便在實際中正確應(yīng)用FFT。二、實驗步驟1．復(fù)習(xí)DFT的定義、性質(zhì)和用DFT作譜分析的有關(guān)內(nèi)容。復(fù)習(xí)FFT算法原理與編程思想，并對照DIT—FFT運算流圖和程序框圖，2．讀懂本實驗提供的FFT子程序。3．編制信號產(chǎn)生子程序，產(chǎn)生以下典型信號供譜分析：xnRn 1 4n1，0n3xn8n，4n720，其它n4n，0n3xnn3，4n730，其它nxncosn 4 4xnsinn 5 8xtcos8tcos16tcos20t6應(yīng)當(dāng)注意，如果給出的是連續(xù)信號xt，則首先要根據(jù)其最高頻率確定采樣速率f以及由 a s頻率分辨率選擇采樣點數(shù)N，然后對其進(jìn)行軟件采樣(即計算xnxnT，a0nN1)，產(chǎn)生對應(yīng)序列xn。對信號xt，頻率分辨率的選擇要以能分辨開其中6的三個頻率對應(yīng)的譜線為準(zhǔn)則。對周期序列，最好截取周期的整數(shù)倍進(jìn)行譜分析，否則有可能產(chǎn)生較大的分析誤差。4．編寫主程序下圖給出了主程序框圖，供參考。本實驗提供FFT子程序和通用繪圖子程序。開始開始讀入長度N調(diào)用信號產(chǎn)生子程序產(chǎn)生實驗信號調(diào)用繪圖子程序(函數(shù))繪制時間序列波形圖調(diào)用FFT子程序(函數(shù))計算信號的DFT調(diào)用繪圖子程序(函數(shù))繪制kX曲線結(jié)束主程序框圖三、實驗結(jié)果直接運行程序，按照實驗內(nèi)容及程序提示鍵入1~8，分別對xn~xn及 1 6xnxnxn、xnxnjxn進(jìn)行譜分析。輸出xn~xn的波形及74584515其8點DFT和16點DFT，xn的16點、32點和64點采樣序列及其DFT。64.實驗程序X1的實驗程序x1=[1111];y11=fft(x1,8);y12=fft(x1,16);subplot(2,3,1);stem(0:3,x1);title('函數(shù)x1的圖像');boxonsubplot(2,2,2);stem(0:7,abs(y11));title('N=8時，x1的DFT圖像');boxonsubplot(2,2,4);stem(0:15,abs(y12));title('N=16時，x1的DFT圖像');boxonX2的實驗程序xa=1:1:4;xb=4:-1:1;x2=[xaxb];y21=fft(x2,8);y22=fft(x2,16);subplot(2,2,1);stem(0:7,x2);title('函數(shù)x2的圖像');subplot(2,2,2);boxonstem(0:7,abs(y21));title('N=8時，x2的DFT圖像');boxonsubplot(2,2,4);stem(0:15,abs(y22));title('N=16時，x1的DFT圖像');boxonX3的實驗程序xa=4:-1:1;xb=1:1:4;x3=[xaxb];y31=fft(x2,8);y32=fft(x2,16);subplot(2,2,1);stem(0:7,x3);title('函數(shù)x3的圖像');subplot(2,2,2);stem(0:7,abs(y31));title('N=8時，函數(shù)x3的DFT圖像');subplot(2,2,3);stem(0:15,abs(y32));title('N=16時，函數(shù)x3的DFT圖像');title('N=16時，x2的DFT圖像');boxonX4X5X7X8的實驗程序subplot(441)n=0:1:7x4=cos(0.25*pi*n)stem(n,x4)xlabel('n')ylabel('x4(n)')title('x4(n)')subplot(442)n=0:1:15x44=cos(0.25*pi*n)stem(n,x44)xlabel('n')ylabel('x44(n)')title('x44(n)')subplot(443)n=0:1:7stem(n,abs(fft(x4)))xlabel('n')ylabel('|X4(k)|')title('8pointFFT')subplot(444)n=0:1:15stem(n,abs(fft(x44)))xlabel('n')ylabel('|X44(k)|')title('16pointFFT')subplot(445)n=0:1:7x5=sin(0.125*pi*n)stem(n,x5)xlabel('n')ylabel('x5(n)')title('x5(n)')subplot(446)n=0:1:15x55=sin(0.125*pi*n)stem(n,x55)xlabel('n')ylabel('x55(n)')title('x55(n)')subplot(447)n=0:1:7stem(n,abs(fft(x5)))xlabel('n')ylabel('|X5(k)|')title('8pointFFT')subplot(448)n=0:1:15stem(n,abs(fft(x55)))subplot(449)n=0:1:7x7=cos(0.25*pi*n)+sin(0.125*pi*n)stem(n,x7)xlabel('n')ylabel('x7(n)')title('x7(n)')subplot(4,4,10)n=0:1:15x77=cos(0.25*pi*n)+sin(0.125*pi*n)stem(n,x77)xlabel('n')ylabel('x77(n)')title('x77(n)')subplot(4,4,11)n=0:1:7stem(n,abs(fft(x7)))xlabel('n')ylabel('|X7(k)|')title('8pointFFT')subplot(4,4,12)n=0:1:15stem(n,abs(fft(x77)))xlabel('n')ylabel('|X77(k)|')title('16pointFFT')subplot(4,4,13)n=0:1:7x8=cos(0.25*pi*n)+j*sin(0.125*pi*n)stem(n,x8)xlabel('n')ylabel('x8(n)')title('x8(n)')subplot(4,4,14)n=0:1:15x88=cos(0.25*pi*n)+j*sin(0.125*pi*n)stem(n,x88)xlabel('n')ylabel('x88(n)')title('x88(n)')subplot(4,4,15)n=0:1:7stem(n,abs(fft(x8)))xlabel('n')ylabel('|X8(k)|')title('8pointFFT')subplot(4,4,16)n=0:1:15stem(n,abs(fft(x88)))xlabel('n')ylabel('|X88(k)|')title('16pointFFT')X6的實驗程序subplot(231)n=0:16fs=64x6=cos(8*pi*n/fs)+cos(16*pi*n/fs)+cos(20*pi*n/fs)stem(n,x6)xlabel('n')ylabel('x6(n)')title('x6(n)')subplot(232)n=0:32fs=64x66=cos(8*pi*n/fs)+cos(16*pi*n/fs)+cos(20*pi*n/fs)stem(n,x66)holdonxlabel('n')ylabel('x66(n)')title('x66(n)')subplot(233)n=0:64fs=64x666=cos(8*pi*n/fs)+cos(16*pi*n/fs)+cos(20*pi*n/fs)stem(n,x666)holdonxlabel('n')ylabel('x666(n)')title('x666(n)')subplot(234)stem(abs(fft(x6)))holdonxlabel('n')ylabel('|X6(k)|')title('16pointFFT')subplot(235)stem(abs(fft(x66)))holdonxlabel('n')ylabel('|X66(k)|')title('32pointFFT')subplot(236)stem(abs(fft(x666)))holdonxlabel('n')ylabel('|X666(k)|')title('64pointFFT')1．xn及其8點和16點DFT122．nx2及其8點和16點DFT3．3．nx3及其8點和16點DFT4．nx4的8點和16點波形及其DFT5．5．nx5的8點和16點波形及其DFT6．nx6的16點、32點和64點采樣序列波形及其DFT77．①計算并圖示nRnxnxnx8547][和nRnxnxnx16547][及其DFT。②程序自動計算并繪圖驗證DFT的共軛對稱性。當(dāng)N=16時，nNxnx44，nNxnx55。即nx4為nx7的共軛對稱分量，而nx5是nx7的共軛反對稱分量。根據(jù)DFT的共軛對稱性，應(yīng)有以下結(jié)果：]Im[]Re[][7777kXjkXnxDFTkXxnxnxn的8點和16點波形及其DFT7 4 5opxn的8點和16點DFT8程序計算結(jié)果如下：程序計算結(jié)果如下：kNXkXkXe*88821及][884kXIDFTnxnxer，正好與圖中xn的16點Xk及xn相同。4 4 4jXk1XkX*Nk及xnxnIDFT[Xk]/j，正好與圖中 8o 28 85 8i 8o16點的Xk及xn相同。55四、實驗總結(jié)本實驗主要是求xn、xn、xn、xn、xn、xn、xn、xn的DFT 1 2 3 4 5 6 7 8變換。其中xn~xn是直接給出了離散序列，而xn、xn則是經(jīng)過xn、xn 1 6 7 8 4 5運算得到的xnxnxn、xnxnjxn。離散傅立葉變換可以看作是 7 4 5 8 4 5Xz在zej2Nk時的Z變換，即表明xn的N點DFT是xn的Z變換在單位圓上的N點等間隔采樣。離散傅立葉變換也可以看作Xk在w2k時的傅立葉變換，即表明XkN可以看作xn的傅立葉變換Xejw在區(qū)間[0,2]上的N點等間隔采樣。五、思考題1．在N=8時，xn和xn的幅頻特性會相同嗎？為什么？N=16呢？ 2 3答：N=8時兩個的幅頻特性相同，因為其不為0的區(qū)間長度正好是8。N=16時兩個的幅頻特性不相同。2、如果周期信號的周期預(yù)先不知道，如何用FFT進(jìn)行譜分析？答：可以先求出它的離散傅立葉變換，通過其DFT的圖形來確定出原始序列的周期。實驗二：用窗函數(shù)法設(shè)計FIR數(shù)字濾波器實驗?zāi)康模?1)熟悉矩形窗、海寧窗、漢明窗和布萊克曼窗。掌握用上述窗函數(shù)法設(shè)計FIR數(shù)字濾波器的原理和方法。熟悉線性相位FIR數(shù)字濾波器特性。了解各種窗函數(shù)對濾波特性的影響。實驗原理和方法：如果所希望的濾波器的理想頻率響應(yīng)函數(shù)為Hd(ejω)，則其對應(yīng)的單位脈沖響應(yīng)為1H(ej)ejndh(n)用窗函數(shù)w(n)將h(n)截斷并進(jìn)處，得到:dh(n)h(n)(n)dh(n)就作為實際設(shè)計的FIR數(shù)字濾波器的單位脈沖響應(yīng)序列，其頻率響應(yīng)函數(shù)H(ejω)為H(ej)N1h(n)ejnn0如果要求線性相位特性，則h(n)還必須滿足：h(n)h(N1n)根據(jù)上式中的正、負(fù)號和長度N的奇偶性又將線性相位FIR濾波器分成四類。要根據(jù)所設(shè)計的濾波特性正確選擇其中一類。例如，要設(shè)計線性相位低通特性，可選擇h(n)=h(N-1-n)一類，而不能選h(n)=-h(N-1-n)一類。實驗內(nèi)容及步驟：復(fù)習(xí)用窗函數(shù)法設(shè)計FIR數(shù)字濾波器一節(jié)內(nèi)容，閱讀本實驗原理，掌握設(shè)計步驟。編寫程序。編寫能產(chǎn)生四種窗函數(shù)的子程序。編寫主程序。實驗程序及結(jié)果：（1）矩形窗clearall;n=0:1:14;wR=ones(1,15);%編寫矩形窗hd=sin(0.25*pi*(n-7+eps))./(pi*(n-7+eps));%讀入hd（n）函數(shù)h1=hd.*wR;%計算h（n）N=64;H1=fft(h1,N);%調(diào)用子程序計算H（k）n=0:N-1;w=2*pi/64*n;subplot(2,2,1);plot(w,fftshift(20*log10((abs(H1)))));%畫幅度曲線gridxlabel('w/rad')ylabel('20lg|H(jw)|/dB')title('幅度曲線和相頻曲線（n=15）');n=0:N-1;w=2*pi/64*n;subplot(2,2,3);plot(w,unwrap(phase(H1)));%畫相頻曲線gridxlabel('w/rad')clearall;n=0:1:32;wR=ones(1,33);%編寫矩形窗hd=sin(0.25*pi*(n-16+eps))./(pi*(n-16+eps));%讀入hd（n）函數(shù)h1=hd.*wR;%計算h（n）N=64;H1=fft(h1,N);%調(diào)用子程序計算H（k）n=0:N-1;w=2*pi/64*n;subplot(2,2,2);plot(w,fftshift(20*log10((abs(H1)))));%畫幅度曲線gridxlabel('w/rad')ylabel('20lg|H(jw)|/dB')title('幅度曲線和相頻曲線（n=33）');n=0:N-1;w=2*pi/64*n;subplot(2,2,4);plot(w,unwrap(phase(H1)));%畫相頻曲線gridxlabel('w/rad')結(jié)果圖像：（2）漢寧窗clearall;n=0:1:14;wH=0.5*(1-cos(2*pi/14*n));%編寫漢寧窗hd=sin(0.25*pi*(n-7+eps))./(pi*(n-7+eps));%讀入hd（n）函數(shù)h1=hd.*wH;%計算h（n）N=64;H1=fft(h1,N);%調(diào)用子程序計算H（k）n=0:N-1;w=2*pi/64*n;subplot(2,2,1);subplot(2,2,1);plot(w,fftshift(20*log10((abs(H1)))));%畫幅度曲線gridxlabel('w/rad')ylabel('20lg|H(jw)|/dB');title('幅度曲線和相頻曲線（n=15）');n=0:N-1;w=2*pi/64*n;subplot(2,2,1);subplot(2,2,3);plot(w,unwrap(phase(H1)));%畫相頻曲線gridxlabel('w/rad')n=0:1:32;wH=0.5*(1-cos(2*pi/32*n));%編寫漢寧窗hd=sin(0.25*pi*(n-16+eps))./(pi*(n-16+eps));%讀入hd（n）函數(shù)h1=hd.*wH;%計算h（n）N=64;H1=fft(h1,N);%調(diào)用子程序計算H（k）n=0:N-1;w=2*pi/64*n;subplot(2,2,1);subplot(2,2,2);plot(w,fftshift(20*log10((abs(H1)))));%畫幅度曲線gridxlabel('w/rad')ylabel('20lg|H(jw)|/dB')title('幅度曲線和相頻曲線（n=33）');n=0:N-1;w=2*pi/64*n;subplot(2,2,1);subplot(2,2,4);plot(w,unwrap(phase(H1)));%畫相頻曲線gridxlabel('w/rad')結(jié)果圖像：（3）海明窗：clearall;n=0:1:14;wH=0.54-0.46*cos(2*pi*n/(14+eps));%編寫海明窗hd=sin(0.25*pi*(n-7+eps))./(pi*(n-7+eps));%讀入hd（n）函數(shù)h1=hd.*wH;%計算h（n）N=64;H1=fft(h1,N);%調(diào)用子程序計算H（k）n=0:N-1;w=2*pi/64*n;subplot(2,2,1);subplot(2,2,1);plot(w,fftshift(20*log10((abs(H1)))));%畫幅度曲線gridxlabel('w/rad')ylabel('20lg|H(jw)|/dB')title('幅度曲線和相頻曲線（n=15）');n=0:N-1;w=2*pi/64*n;subplot(2,2,1);subplot(2,2,3);plot(w,unwrap(phase(H1)));%畫相頻曲線gridxlabel('w/rad')n=0:1:32;wH=0.54-0.46*cos(2*pi*n/(32+eps));%編寫海明窗hd=sin(0.25*pi*(n-16+eps))./(pi*(n-16+eps));%讀入hd（n）函數(shù)h1=hd.*wH;%計算h（n）N=64;H1=fft(h1,N);%調(diào)用子程序計算H（k）n=0:N-1;w=2*pi/64*n;subplot(2,2,1);subplot(2,2,2);plot(w,fftshift(20*log10((abs(H1)))));%畫幅度曲線gridxlabel('w/rad')ylabel('20lg|H(jw)|/dB')title('幅度曲線和相頻曲線（n=33）');n=0:N-1;w=2*pi/64*n;subplot(2,2,1);subplot(2,2,4);plot(w,unwrap(phase(H1)));%畫相頻曲線gridxlabel('w/rad')結(jié)果圖像(4)布萊克曼窗n=0:1:14;wB=0.42-0.5*cos(2*pi/(14+eps)*n)+0.08*cos(4*pi/(14+eps)*n);%編寫布萊克曼窗hd=sin(0.25*pi*(n-7+eps))./(pi*(n-7+eps));%讀入hd（n）函數(shù)h1=hd.*wB;%計算h（n）N=64;H1=fft(h1,N);%調(diào)用子程序計算H（k）n=0:N-1;w=2*pi/64*n;subplot(2,2,1);subplot(2,2,1);plot(w,fftshift(20*log10((abs(H1)))));%畫幅度曲線gridxlabel('w/rad')ylabel('20lg|H(jw)|/dB')title('幅度曲線和相頻曲線（n=15）');n=0:N-1;w=2*pi/64*n;subplot(2,2,1);subplot(2,2,3);plot(w,unwrap(phase(H1)));%畫相頻曲線gridxlabel('w/rad')n=0:1:32;wB=0.42-0.5*cos(2*pi/(32+eps)*n)+0.08*cos(4*pi/(32+eps)*n);%編寫布萊克曼窗hd=sin(0.25*pi*(n-16+eps))./(pi*(n-16+eps));%讀入hd（n）函數(shù)h1=hd.*wB;%計算h（n）N=64;H1=fft(h1,N);%調(diào)用子程序計算H（k）n=0:N-1;w=2*pi/64*n;subplot(2,2,1);subplot(2,2,2);plot(w,fftshift(20*log10((abs(H1)))));%畫幅度曲線gridxlabel('w/rad')ylabel('20lg|H(jw)|/dB')title('幅度曲線和相頻曲線（n=33）');n=0:N-1;w=2*pi/64*n;subplot(2,2,1);subplot(2,2,4);plot(w,unwrap(phase(H1)));%畫相頻曲線gridxlabel('w/rad')結(jié)果圖像：實驗分析：窗函數(shù)法特點：窗口法設(shè)計的主要優(yōu)點是簡單，使用方便。窗口函數(shù)大多有封閉的公式可循，性能、參數(shù)都已有表格、資料可供參考，計算程序簡便，所以很實用。缺點是通帶和阻帶的截止頻率不易控制。窗口長度N和窗函數(shù)類型對濾波特性的影響：增加窗的長度，可以減少窗的主瓣寬度，從而減少H()過渡帶的帶寬，但增加N并不能減少帶內(nèi)波動以及加大阻帶。這兩個指標(biāo)只能從窗函數(shù)的形狀上找解決方法。常用的6種窗口的基本參數(shù)及性能窗類型窗類型窗譜性能指標(biāo)加窗后濾波器的性能指標(biāo)最大旁瓣(dB)主瓣帶寬阻帶衰減dB)(過渡帶寬Δω矩形窗-132*2π/N-210.9*2π/N三角窗-254*2π/N-252.1*2π/Nhanning窗-314*2π/N-443.1*2π/Nhamming窗-414*2π/N-533.3*2π/NBlackman窗-576*2π/N-745.5*2π/NKaiser窗（β=7.865）-57.07-805.0*2π/N通過對窗函數(shù)選擇，可以獲得符合要求的阻帶衰減指標(biāo)。實驗總結(jié)及心得用窗函數(shù)法設(shè)計FIR濾波器的主要特點：設(shè)H(ejw)＝FT[h(n)]為希望逼近的頻響特性函數(shù),H(ejw)=FT[h(n)]為用窗函數(shù)法 d d設(shè)計的實際濾波器的頻響函數(shù)。通常取H(ejw)相應(yīng)的理想頻響特性作為H(ejw)。d知識要點如下：1．希望逼近的理想濾波器頻響函數(shù)H(ejw)的表達(dá)式。因為數(shù)字濾波器一般要求設(shè)計d成線性相位特性，所以H(ejw)必須滿足上述線性相位FIR濾波器的頻域