數(shù)字信號處理實習(xí)報告

1、 中國地質(zhì)大學(xué)（武漢） 數(shù)字信號處理上機實習(xí) 學(xué)生姓名： 班 級：071132 學(xué) 號：2013100 指導(dǎo)老師：王曉莉 題目一 離散卷積計算1、 實驗題目設(shè)線性時不變（LTI）系統(tǒng)的沖激響應(yīng)為h(n)，輸入序列為x(n)1、h(n)=(0.8)n，0n4； x(n)=u(n)-u(n-4)2、h(n)=(0.8)nu(n), x(n)=u(n)-u(n-4)3、h(n)=(0.8)nu(n), x(n)=u(n)求以上三種情況下系統(tǒng)的輸出y(n)，顯示輸入和輸出波形。2、 實驗?zāi)康?.理解和掌握離散卷積計算；2.學(xué)習(xí)如何用Mtalab實現(xiàn)離散卷積計算。三、算法設(shè)計離散卷積定義為：1、 ， (

2、a) 當(dāng) 時，；(b) 當(dāng)時,(0.8)n ；(c) 當(dāng)時，(0.8)n；(d) 當(dāng)時，；2、， (a) 當(dāng) 時，；(b) 當(dāng)時,(0.8)n；(c) 當(dāng)時，(0.8)n；(d) 當(dāng)時，(0.8)n；(e) 當(dāng)時，；3、， (a) 當(dāng) 時，；(b) 當(dāng)時,(0.8)n；(c) 當(dāng)時，(0.8)n；(d) 當(dāng)時，；4、 程序分析所用到的函數(shù)：（1）y=conv（x.,h）：卷積運算函數(shù)，計算;（2）n1=0:4：n1取04；（3）subplot(m,n,p):subplot()函數(shù)是將多個圖畫到一個平面上的工具。其中，m表示是圖排成m行，n表示圖排成n列，也就是整個figure中有n個圖是排成一

3、行的，一共m行，如果m=2就是表示2行圖。p表示圖所在的位置，p=1表示從左到右從上到下的第一個位置。（4）title(content ):title()函數(shù)的功能是為當(dāng)前坐標(biāo)系添加標(biāo)題“content ”。五、程序設(shè)計n=0:4;h=0.8.n;x=1 1 1 1subplot(331);stem(x); title('x(n)');y=conv(x,h);subplot(332);stem(h); title('h(n)');subplot(333);stem(y); title('y(n)'); n=0:40;h=0.8.n;x=1 1 1

4、 1subplot(334);stem(x); title('x(n)');y=conv(x,h);subplot(335);stem(h); title('h(n)');subplot(336);stem(y); title('y(n)'); n=0:40;h=0.8.n;x=zeros(1,0),ones(1,40);subplot(337);stem(x); title('x(n)');y=conv(x,h);subplot(338);stem(h); title('h(n)');subplot(339);s

5、tem(y); title('y(n)');六、運行結(jié)果圖中從左至右三列依次對應(yīng)、及卷積結(jié)果題目二 離散傅立葉變換1、 實驗題目設(shè)有離散序列 分析下列三種情況下的幅頻特性。（1） 采集數(shù)據(jù)長度N=16，分析16點的頻譜，并畫出幅頻特性。采集數(shù)據(jù)長度N=16，并補零到64點，分析其頻譜，并畫出幅頻特性。（2） 采集數(shù)據(jù)長度N=64，分析64點的頻譜，并畫出幅頻特性。觀察三幅不同的幅頻特性圖，分析和比較它們的特點及形成原因。2、 實驗?zāi)康?、理解掌握DFT及FFT算法；2、利用FFT算法計算信號的頻譜。三、算法設(shè)計當(dāng)抽樣數(shù)N=2M時，以下為蝶形算法圖。1、 當(dāng)N=2M時，則要進(jìn)行M

6、次分解，即進(jìn)行M級蝶形單元的計算；2、按自然順序輸入，輸出是碼位倒置；3、每一級包含N/2個基本蝶形運算；4、第L級有2L-1個蝶群，蝶群間隔為N/2L-1；四、程序分析所用到的函數(shù)：（1） fft（）：求x的一維傅里葉變換，計算X（k）；使用方法Xk=fft(xn)，例如：N=8;n=0:N-1;xn=4 3 2 6 7 8 9 0;Xk=fft(xn)（2）stem（）:顯示函數(shù)圖像。（3）xlabel():橫坐標(biāo)軸的名稱；（4）ylabel():縱坐標(biāo)軸的名稱；五、程序設(shè)計%x(n)=cos(0.48n)+cos(0.52n)n=0:15; %產(chǎn)生序列x(n) 取16點x=cos(0.4

7、8*pi*n)+cos(0.52*pi*n);subplot(231);stem(n,x);title('采集數(shù)據(jù)長度N=16'); %顯示x(n)x1k=fft(x,16);%進(jìn)行16點傅里葉變換subplot(234);stem(0:15,abs(x1k),'.'); %顯示X（k）xlabel('n');ylabel('X1(k)'); title('16點傅立葉變換頻譜'); x=x zeros(1,48); %取16點，補零到64點x2k=fft(x,64);%進(jìn)行64點傅里葉變換subplot(232)

8、;stem(x);title('采集N=16,補零到64');subplot(235);stem(0:63,abs(x2k),'.');xlabel('n');ylabel('X2(k)');title('64點傅立葉變換頻譜'); n=0:63; %產(chǎn)生序列x(n)，取64點x=cos(0.48*pi*n)+cos(0.52*pi*n);subplot(233);stem(n,x);title('采集數(shù)據(jù)長度N=64');x3k=fft(x);%做64點傅里葉變換subplot(236);stem

9、(0:63,abs(x3k),'.');xlabel('n');ylabel('X3(k)');title('64點傅立葉變換頻譜');六、運行結(jié)果七結(jié)果分析N點DFT的頻譜分辨率是2 /N。一節(jié)指出可以通過補零觀察到更多的頻點，但是這并不意味著補零能夠提高真正的頻譜分辨率。這是因為xn 實際上是x(t) 采樣的主值序列，而將xn補零得到的x'n 周期延拓之后與原來的序列并不相同，也不是x(t) 的采樣。因此是不同離散信號的頻譜。對于補零至M點的x'的DFT，只能說它的分辨率2 /M僅具有計算上的意義，并不是真正的

10、、物理意義上的頻譜。頻譜分辨率的提高只能通過提高采樣頻率實現(xiàn)。第三幅圖形取樣點數(shù)較多，分辨率最高，而第一幅圖形則顯得較為稀疏，主要是因為取樣點數(shù)太少。題目三 IIR濾波器的設(shè)計1、 實驗題目1、設(shè)計一個切比雪夫數(shù)字低通濾波器，設(shè)計指標(biāo)如下： 通帶截止頻率：0.2，幅度衰減不大于3分貝 阻帶截止頻率：0.3，幅度衰減大于20分貝2、分析不同濾波器的特點和結(jié)果。3、編程設(shè)計實現(xiàn)IIR濾波器。二、實驗?zāi)康?.理解和掌握不同IIR濾波器的性質(zhì)、特點。2.掌握用Matlab設(shè)計切比雪夫濾波器和橢圓濾波器的方法。三、算法設(shè)計1、濾波器類型切比雪夫濾波器：切比雪夫濾波器，又名“車比雪夫濾波器”，是在通帶或阻

11、帶上頻率響應(yīng)幅度等波紋波動的濾波器。2、變換方法（a）沖激響應(yīng)不變法沖激響應(yīng)不變法的基本原理是從濾波器的沖激響應(yīng)出發(fā)，對模擬濾波器沖激響應(yīng)h(t)進(jìn)行取樣，所得到的離散序列h(nT)作為數(shù)字濾波器的單位取樣響應(yīng)。H(z)是由H(s)通過下式的對應(yīng)關(guān)系得到。（b）雙線性變換是在所得到滿足性能指標(biāo)要求的模擬濾波器的基礎(chǔ)上,通過變換,從而得到相應(yīng)的數(shù)字濾波器。4、 程序分析所用到的函數(shù)：（1）N,wn=buttord(Wp,Ws,ap,as,'s')：階數(shù)，截止頻率=buttord（通帶截止頻率，阻帶截止頻率，通帶衰減，阻帶衰減，模擬）；（2）num,den=bilinear(b,a

12、,fs)：雙線性變換可選頻率預(yù)畸變，在IIR數(shù)字濾波器設(shè)計中有把模擬的原型濾波器轉(zhuǎn)變成相應(yīng)的數(shù)字濾波器，即通過已知S變換中模擬濾波器的系數(shù)求數(shù)字濾波器的系數(shù)。其中有兩種方法來進(jìn)行這種變換，一種是激沖不變法，另一種是雙線性Z變換法。bilinear 函數(shù)便是雙線性Z變換法。函數(shù)中的b和a是原模擬濾波器S變換中的分子和分母系數(shù)，num和den是數(shù)字濾波器的系數(shù)。（3）grid on：增加主要網(wǎng)格線為當(dāng)前軸；（4）h,w=freqz(num,den)：求頻率響應(yīng)。函數(shù)的輸出：a.濾波器的頻率響應(yīng)H（N點） b.頻率向量W(N點，且單位為弧度)，頻率向量W是均勻分布在濾波器的上半?yún)^(qū)，即：0：

13、pi，這些點上的頻率響應(yīng)都將通過此函數(shù)計算出來。；濾波器的系數(shù)：分子為num，分母為den。    五、程序設(shè)計clc;Rs=20; Rp=3;Wp1=0.2*pi;Ws1=0.3*pi;N,Wn=cheb1ord(Wp1,Ws1,Rp,Rs,'s');num1,den1=cheby1(N,Rp,Wn,'s');num,den=bilinear(num1,den1,0.5); subplot(2,2,1);zplane(num,den);title('零極點圖')w=0:pi/256:pi;h=freqz(num,d

14、en,w);g=abs(h);g1=angle(h);subplot(223);plot(w/pi,g);gridaxis(0 1 0 1);xlabel('頻率'); ylabel('幅度'); title('幅頻響應(yīng)');subplot(224);plot(w/pi,g1);grid onaxis(0 1 -5 5);xlabel('頻率'); ylabel('相位'); title('相頻響應(yīng)');%測試t=1:300;I=sin(0.1*pi*t)+sin(0.8*pi*t)+sin(0.9

15、5*pi*t);%設(shè)計正弦波figure;subplot(311);plot(I);title('原信號x(n)=sin(0.1*pi*t)+sin(0.8*pi*t)+sin(0.95*pi*t)');A=filter(num,den,I);%正弦波通過濾波器subplot(312);plot(A);title('濾波后y(n)');6、 運行結(jié)果（1） 零極點圖、幅頻響應(yīng)、相頻響應(yīng)（2）輸入為I=sin(0.1*pi*t)+sin(0.8*pi*t)+sin(0.95*pi*t)時的濾波結(jié)果7、 結(jié)果分析由幅頻特性可知，橢圓濾波器在通帶內(nèi)的頻率特性是平坦的,

16、并且隨著頻率的增加而衰減。通帶截止頻率：0.2，幅度衰減不大于3分貝，阻帶截止頻率：0.3，幅度衰減大于20分貝。正弦信號在經(jīng)過IIR濾波器濾波后，由傅里葉變換后的頻譜圖可看出高頻信號被濾除，低頻信號被保留了下來。題目四 FIR濾波器的設(shè)計一、實驗內(nèi)容選取合適窗函數(shù)設(shè)計一個線性相位FIR低通濾波器，使它滿足如下性能指標(biāo)：通帶截止頻率：p=0.66,通帶截止頻率處的衰減不大于3分貝；阻帶截止頻率：s=0.5，阻帶衰減不小于40分貝二、實驗?zāi)康?、掌握用窗函數(shù)法設(shè)計FIR濾波器的原理和方法。2、熟悉線性相位濾波器特性。3、了解各種窗函數(shù)對濾波器特性的影響。三、原理及算法概要1、算法通過其通帶截止頻

17、率p與阻帶截止頻率s算出其過渡帶的寬度與濾波器的長度，從而得到理想濾波器的截止頻率，根據(jù)所要求的理想濾波器，得到hd（n）。由于其通帶截止頻率處的衰減不大于3分貝與阻帶衰減不小于40分貝，我選擇最接近的漢寧窗，最后調(diào)用函數(shù)h=hd.*win 及freqz(h,1,512)得到實際漢寧窗的響應(yīng)和實際濾波器的幅度響應(yīng)。2、原理概要利用窗函數(shù)法設(shè)計FIR 濾波器，F(xiàn)IR濾波器的最大特點是其相位特性可以設(shè)計為嚴(yán)格的線性，而其幅值可以任意設(shè)置，這樣輸出波形就不會相位失真。理想低通濾波器的單位取樣響應(yīng)hd（n）是無限長的，所以要用一個有限長的因果序列h（n）進(jìn)行逼近，最有效的方法是截斷hd（n），即用有限

18、長的窗函數(shù)w（n）來截取hd（n），表示為h（n）=hd（n）w（n）。為獲得線性相位的FIR濾波器，h（n）必須滿足中心對稱條件，序列h（n）應(yīng)有一定的延遲，且=（N-1）/2頻率響應(yīng)逼近hd（ejw）的FIR濾波器，最簡單的窗函數(shù)為矩形窗： 1 n<(N-1)/2 W（n）=0 n>(N-1)/2加窗后的頻譜 加窗后使實際頻響偏離理想頻響，影響主要有兩個方面：（1） 通帶和阻帶之間存在過渡帶，過渡帶寬度取決于窗函數(shù)頻響的主瓣寬度。（2） 通帶和阻帶區(qū)間有紋波，這是由窗函數(shù)的旁瓣引起的，旁瓣越多，紋波越多。增加窗函數(shù)的寬度N，其主瓣寬度減小，但不改變旁瓣的相對值。為了改善濾波器的

19、性能，要求窗函數(shù)的主瓣寬度盡可能窄，以獲得較窄的過渡帶；旁瓣衰減盡可能大，數(shù)量盡可能大，從而改善紋波狀況，使實際頻響H（ej）更好地逼近理想頻響Hd（ej）。除了矩形窗外，一般還可以采用以下幾種窗函數(shù) 漢寧窗：海明窗布來克曼窗4、 程序設(shè)計 clc; wp=0.66*pi;ws=0.5*pi;wd=wp-ws;%hanning窗N_hann=ceil(8*pi/wd); wc=(wp+ws)/2;if mod(N_hann,2)=0 N_hann=N_hann+1;endhann_win=hanning(N_hann);b_hann=fir1(N_hann-1,wc/pi,'high&

20、#39;,hann_win);freq_axis=0:pi/512:pi-pi/512;freqn=0:511/512;H_hann=freqz(b_hann,1,512);figure(1);subplot(311);plot(freqn,20*log10(abs(H_hann);xlabel('歸一化頻率w/pi');ylabel('幅度譜（dB）');title('采用hanning窗設(shè)計fir濾波器幅度譜');subplot(312);plot(freqn,angle(H_hann);xlabel('歸一化頻率w/pi')

21、;ylabel('相位譜');title('采用hanning窗設(shè)計fir濾波器相位譜');subplot(313);stem(b_hann);xlabel('n');ylabel('h(n)');title('單位脈沖響應(yīng)序列');%hamming窗N_hamm=ceil(8*pi/wd);wc=(wp+ws)/2;if mod(N_hamm,2)=0 N_hamm=N_hamm+1;endhamm_win=hamming(N_hamm);b_hamm=fir1(N_hamm-1,wc/pi,'high&

22、#39;,hamm_win);H_hamm=freqz(b_hamm,1,512);figure(2); subplot(311);plot(freqn,20*log10(abs(H_hamm);xlabel('歸一化頻率w/pi');ylabel('幅度譜（dB）');title('采用hamming窗設(shè)計fir濾波器幅度譜');subplot(312);plot(freqn,angle(H_hamm);xlabel('歸一化頻率w/pi');ylabel('相位譜');title('采用hamming窗

23、設(shè)計fir濾波器相位譜');subplot(313);stem(b_hann);xlabel('n');ylabel('h(n)');title('單位脈沖響應(yīng)序列');5 程序運行結(jié)果（1）采用hanning窗設(shè)計FIR濾波器（2）采用hamming窗設(shè)計FIR濾波器6 結(jié)果分析上圖為選用漢寧窗和漢明窗的FIR 濾波器的相頻特性、幅頻特性及單位取樣響應(yīng)hd（n）。由運行得出的幅頻特性可知，滿足設(shè)計指標(biāo)通帶截止頻率：p=0.66；阻帶截止頻率：s=0.5。得出的窗主瓣寬度較寬，旁瓣寬度較窄，也合乎要求。題目五一、實驗內(nèi)容錄制一段電話的撥號

24、音，并對錄制的信號進(jìn)行采樣；畫出采樣后信號的時域波形和頻譜圖；必要的話對采集的信號進(jìn)行濾波，畫出濾波后信號的時域波形和頻譜，并對濾波前后的信號進(jìn)行對比，根據(jù)信號的頻譜解析出撥號的數(shù)字。最后，用MATLAB設(shè)計一信號處理系統(tǒng)界面。二、實驗?zāi)康?、掌握用濾波器處理語音信號的原理和方法。2、掌握用matlab設(shè)計界面。三、原理及算法概要調(diào)用函數(shù)function pushbutton1_Callback(hObject, eventdata, handles)實現(xiàn)一個信號處理系統(tǒng)界面。選擇button1時，用雙線性變換法設(shè)計濾波器來對信號進(jìn)行處理。讀取語音信號，對語音信號進(jìn)行f=8000的頻率進(jìn)行采樣

25、，調(diào)用函數(shù)y1=fft(x1,2048)對所采集的點做2048點FFT變換。先設(shè)計butterworth模擬濾波器，再用雙線性變換法實現(xiàn)模擬濾波器到數(shù)字濾波器的轉(zhuǎn)換。最后調(diào)用函數(shù)f1=filter(bz,az,x2)對加了噪聲的語音信號進(jìn)行濾波，得到濾波后的頻譜圖。4 程序% - Executes on button press in pushbutton1.function pushbutton1_Callback(hObject, eventdata, handles)% hObject handle to pushbutton1 (see GCBO)% eventdata reserve

26、d - to be defined in a future version of MATLAB% handles structure with handles and user data (see GUIDATA)clc;fs=8000; %語音信號采樣頻率為8000x1=wavread('D:matlab2011workshua.wav');t=(0:length(x1)-1)/8000;y1=fft(x1,2048); %對信號做2048點FFT變換f=fs*(0:1023)/2048;plot(handles.axes1,t,x1); %做原始信號的時域波形grid on

27、;axis tight;title(handles.axes1,'原始語音信號');xlabel(handles.axes1,'time(s)');ylabel(handles.axes1,'幅度');plot(handles.axes2,f,abs(y1(1:1024); %做原始信號的FFT頻譜grid on;axis tight;title(handles.axes2,'原始語音信號的FFT頻譜');xlabel(handles.axes2,'Hz');ylabel(handles.axes2,'幅度

28、');%雙線性變換法設(shè)計的巴特沃斯濾波器A1=0.05;A2=0.10;d=A1*cos(2*pi*3800*t)+A2*sin(2*pi*3600*t)'x2=x1+d;wp=0.8*pi;ws=0.85*pi;Rp=1;Rs=15;Fs=8000;Ts=1/Fs;wp1=2/Ts*tan(wp/2); %將模擬指標(biāo)轉(zhuǎn)換為數(shù)字指標(biāo)ws1=2/Ts*tan(ws/2);N,Wn=buttord(wp1,ws1,Rp,Rs,'s'); %選擇濾波器最小階數(shù)Z,P,K=buttap(N); %創(chuàng)建butterworth模擬濾波器Bap,Aap=zp2tf(Z,P,K);b,a=lp2lp(Bap,Aap,Wn); bz,az=bilinear(b,a,Fs); %用雙線性法實現(xiàn)模擬到數(shù)字的轉(zhuǎn)換H,W=freqz(bz,az); %繪制頻率響應(yīng)曲線plot(handles.axes3,W*Fs/(2*pi),abs(H);grid on