數(shù)字信號處理課程設計--基于 matlab 的音樂信號處理和分析

1、<p>　　《數(shù)字信號處理》課程設計</p><p>　　設計題目：基于 MATLAB 的音樂信號處理和分析</p><p>　　一、課程設計的目的2</p><p>　　二、課程設計基本要求2</p><p>　　三、課程設計內(nèi)容2</p><p>　　1、音樂信號的音譜和頻譜觀察2</p&g

2、t;<p>　　2、音樂信號的抽?。p抽樣）4</p><p>　　3、音樂信號的AM調(diào)制8</p><p>　　4、AM調(diào)制音樂信號的同步解調(diào)11</p><p>　　5、音樂信號的濾波去噪18</p><p>　　6、音樂信號的幅頻濾波和相頻濾波24</p><p><b>　　四

3、、問題討論27</b></p><p>　　1、IIR數(shù)字濾波器和FIR數(shù)字濾波器的比較27</p><p>　　2、音樂信號的音調(diào)與信號的什么特征有關？28</p><p>　　3、音樂信號的音色與信號的什么特征有關？28</p><p>　　4、兩種不同音色的音樂信號疊加混疊后，為何人耳還可以分辨？29</p&

4、gt;<p>　　5、音樂信號的幅度與相位特征對信號有哪些影響？29</p><p><b>　　五、心得體會29</b></p><p><b>　　一、課程設計的目的</b></p><p>　　本課程設計通過對音樂信號的采樣、抽取、調(diào)制、解調(diào)等多種處理過程的理論分析和MATLAB實現(xiàn)，使學生進一步鞏

5、固數(shù)字信號處理的基本概念、理論、分析方法和實現(xiàn)方法；使學生掌握的基本理論和分析方法知識得到進一步擴展；使學生能有效地將理論和實際緊密結合；增強學生軟件編程實現(xiàn)能力和解決實際問題的能力。</p><p>　　二、課程設計基本要求</p><p>　　1 學會 MATLAB 的使用， 掌握MATLAB的基本編程語句。</p><p>　　2 掌握在 Windows 環(huán)境

6、下音樂信號采集的方法。</p><p>　　3 掌握數(shù)字信號處理的基本概念、基本理論和基本方法。</p><p>　　4 掌握 MATLAB 設計 FIR 和 IIR 數(shù)字濾波器的方法。 5 掌握使用MATLAB處理數(shù)字信號、進行頻譜分析、設計數(shù)字濾波器的編程方法。</p><p><b>　　三、課程設計內(nèi)容</b></p>

7、<p>　　1、音樂信號的音譜和頻譜觀察</p><p>　　使用windows下的錄音機錄制一段音樂信號或采用其它軟件截取一段音樂信號（要求：時間不超過5s、文件格式為wav文件）</p><p>　?、?使用wavread語句讀取音樂信號，獲取抽樣率；（注意：讀取的信號是雙聲道信號，即為雙列向量，需要分列處理）；</p><p>　?、?輸出音樂信號的

8、波形和頻譜，觀察現(xiàn)象；</p><p>　?、?使用sound語句播放音樂信號，注意不同抽樣率下的音調(diào)變化，解釋現(xiàn)象。</p><p><b>　　查找help：</b></p><p>　　Wavread格式說明：</p><p>　　[y,fs,b]=wavread(‘語音信號’)，采樣值放在向量y中，fs表示采樣頻

9、率（hz），b表示采樣位數(shù)。</p><p>　　【matlab程序如下】</p><p><b>　　clc</b></p><p><b>　　close all</b></p><p><b>　　clear all</b></p><p>　　[

10、y,fs]=wavread('一生有你');%讀取歌曲</p><p><b>　　size(y)</b></p><p>　　y1=y(:,1); %1聲道</p><p>　　y2=y(:,2); %2聲道</p><p>　　N=length(y);%長度</p><p>

11、;<b>　　n=0:N-1;</b></p><p>　　t=n/fs; %t=nT</p><p>　　w=2*n/N; %2pi在長度N上的平均 是個序列</p><p>　　f1=fft(y1);</p><p>　　f2=fft(y2);%傅里葉變換</p><p>　　fi

12、gure %畫圖像</p><p>　　subplot(2,2,1), plot(t(1:1000),y1(1:1000))</p><p>　　title('1聲道時域圖')</p><p>　　xlabel('t')</p><p>　　ylabel('y1')</p>

13、<p>　　subplot(2,2,3),plot(w,abs(f1)*2/N) </p><p>　　title('1聲道頻域圖')</p><p>　　xlabel('w')</p><p>　　ylabel('Y1(w)')</p><p>　　subplot(2,2,2),p

14、lot(t(1:1000),y2(1:1000))</p><p>　　title('2聲道時域圖')</p><p>　　xlabel('t')</p><p>　　ylabel('y2')</p><p>　　subplot(2,2,4),plot(w,abs(f2)*2/N)</p&

15、gt;<p>　　title('2聲道頻域圖')</p><p>　　xlabel('w')</p><p>　　ylabel('Y2(w)')</p><p>　　wavplay(y1,fs)</p><p>　　wavplay(y1,fs/2) 慢放 wavpla

16、y(y1,fs*2) 快放 </p><p>　　wavylay(y2,fs*2) 快放 wavplay(y2,fs/2) 慢放</p><p>　　【程序運行結果如下圖】：</p><p>　　分析：通過觀察音樂信號的波形和頻譜可知所選取的音樂信號頻譜集中在0~0.5pi之間，同時抽樣頻率為fs=44000</p><p>　

17、　2、音樂信號的抽取（減抽樣）</p><p>　?、?觀察音樂信號頻率上限，選擇適當?shù)某槿￠g隔對信號進行減抽樣（給出兩種抽取間隔，代表混疊與非混疊）；</p><p>　?、?輸出減抽樣音樂信號的波形和頻譜，觀察現(xiàn)象，給出理論解釋； </p><p>　?、?播放減抽樣音樂信號，注意抽樣率的改變，比較不同抽取間隔下的聲音，解釋現(xiàn)象。</p><

18、p><b>　　理論基礎：</b></p><p>　　時域抽樣定理：一個頻譜受限的信號f(t)，如果頻譜只占據(jù)-wm~+wm的圍，</p><p>　　則信號f(t)可以用等間隔的抽樣值唯一的表示。而抽樣間隔必須不大于1/（2*fm）.</p><p>　　頻域抽樣定理：一個頻譜受限的信號f(t)，它集中在-tm~+tm的時間范圍內(nèi)，&

19、lt;/p><p>　　若在頻域中以不大于1/（2*tm）的頻率間隔對f(t)的頻譜F(w)進行抽樣，則抽樣后</p><p>　　頻譜F1(w)可以唯一的表示原信號。</p><p>　　【Matlab程序如下】：</p><p>　　%原信號的頻率上限為0.5pi</p><p>　　clear all;close a

20、ll</p><p>　　[y,fs,bits]=wavread('一生有你');</p><p>　　y1=y(:,1); %取一頻道信號</p><p>　　f1=fft(y1);</p><p>　　N=length(y1);</p><p>　　%減抽樣 $$$$$$$$$$$$$$$$$$

21、$減抽樣使抽樣點數(shù)減少，會使柵欄效應更嚴重</p><p>　　D1=2;D2=16;</p><p><b>　　n=0:N-1;</b></p><p>　　t=n/fs; %t=nT</p><p>　　yd1=y1(1:D1:N); %2倍減抽樣</p><p>　　f

22、yd1=fft(yd1,N);</p><p>　　yd2=y1(1:D2:N);</p><p>　　fyd2=fft(yd2,N); %16倍減抽樣</p><p>　　w=2*n/N;%2pi在長度N上的平均</p><p><b>　　figure</b></p><p>　　subplo

23、t(2,3,1);plot(t(1:1000),y1(1:1000));</p><p>　　title('原信號時域圖')</p><p>　　xlabel('t')</p><p>　　ylabel('y1')</p><p>　　subplot(2,3,4);plot(w,abs(f1))

24、;</p><p>　　title('原信號頻域圖')</p><p>　　xlabel('t')</p><p>　　ylabel('y1')</p><p>　　subplot(2,3,2);plot(t(1:1000),yd1(1:1000));</p><p>　

25、　title('2倍減抽樣后的時域圖')</p><p>　　xlabel('t')</p><p>　　ylabel('yd1')</p><p>　　subplot(2,3,5);plot(w,abs(fyd1));</p><p>　　title('2倍減抽樣后的頻域圖')

26、%出現(xiàn)柵欄效應</p><p>　　xlabel('t')</p><p>　　ylabel('fyd1')</p><p>　　subplot(2,3,3);plot(t(1:1000),yd2(1:1000));</p><p>　　title('10倍減抽樣后的時域圖')</p>

27、;<p>　　xlabel('t')</p><p>　　ylabel('yd2')</p><p>　　subplot(2,3,6);plot(w,abs(fyd2));</p><p>　　title('16倍減抽樣后的頻域圖')%柵欄效應更明顯</p><p>　　xlabe

28、l('t')</p><p>　　ylabel('fyd2')</p><p>　　wavplay(yd1,fs);</p><p>　　wavplay(yd1,fs/D1);</p><p>　　wavplay(yd2,fs);</p><p>　　wavplay(yd2,fs/D2)

29、;</p><p>　　%D倍減抽樣后相當于原抽樣頻率縮小了D倍，故仍按原抽樣頻率播放相當于加快播放</p><p>　　【程序運行結果如下圖】：</p><p>　　分析：通過觀察兩種不同抽樣間隔（2倍頻和16倍頻）下的音樂信號可知，當采用較大的抽樣間隔對音樂信號進行抽樣時，頻譜發(fā)生了混疊，而采用較小的抽樣間隔對音樂信號進行抽樣時，頻譜并未發(fā)生混疊。這是因為，抽樣

30、時頻譜發(fā)生混疊的條件是fs<2fh，即抽樣頻率小于信號頻譜的最高頻率。當采用較大的抽樣間隔時抽樣頻率時fs<2fh，所以發(fā)生混疊，而采用較小的抽樣間隔時抽樣頻率時fs>=2fh，則不會發(fā)生混疊。</p><p>　　當我們播放不同抽樣間隔下的音樂信號時，會發(fā)現(xiàn)大抽樣的音樂信號會伴有雜音并且聲音低沉，而小抽樣的音樂信號和原有的音樂信號幾乎無差別，這間接證明了我們以上理論分析的正確性。</p&

31、gt;<p>　　3、音樂信號的AM調(diào)制</p><p>　?、?觀察音樂信號頻率上限，選擇適當調(diào)制頻率對信號進行調(diào)制（給出高、低兩種調(diào)制頻率）；</p><p>　　② 輸出調(diào)制信號的波形和頻譜，觀察現(xiàn)象，給出理論解釋； </p><p>　?、?播放調(diào)制音樂信號，注意不同調(diào)制頻率下的聲音，解釋現(xiàn)象。</p><p><

32、b>　　理論基礎：</b></p><p>　　信號的調(diào)制過程就是將信號頻譜搬移到任何所需的較高頻率范圍。調(diào)制的實質(zhì)是把各種信號的頻譜搬移，使它們互不重疊的占據(jù)不同的頻率范圍，也即信號分別托付于不同頻率的載波上。具體的調(diào)制原理推導在此不再敘述，僅將結論列出：f(t)=g(t)*cos(w0*t)。由此將信號g(t)的頻譜搬移到（2*n+1）w附近，同時音樂信號的頻譜幅度變?yōu)樵瓉淼?/2。<

33、/p><p>　　如果信號的最高頻譜wh超過了ws/2，則各周期延括分量產(chǎn)生頻譜的交疊，稱為是頻譜的混疊現(xiàn)象。根據(jù)奈奎斯特定律可知，若希望頻譜不會發(fā)生混疊，則fs>=fh。</p><p>　　【Matlab程序如下】：</p><p>　　%原信號的頻率上限為0.5pi</p><p><b>　　clc</b>&l

34、t;/p><p><b>　　clear all</b></p><p>　　[y,fs]=wavread('一生有你');</p><p>　　y1=y(:,1);</p><p>　　N=length(y);</p><p>　　n=0:(N-1);</p><

35、p><b>　　%低調(diào)制頻率</b></p><p>　　y2=cos(n*0.5*pi);%調(diào)制頻率為0.5*pi；</p><p>　　y3=y1.*y2';%調(diào)制后的信號；</p><p>　　f1=fft(y1);%對原信號做fft變換；</p><p>　　f3=fft(y3);%對調(diào)制后的信號做

36、fft變換；</p><p>　　t=n/fs;%t=nT</p><p>　　N=length(y);</p><p>　　w=2*n/length(y);</p><p><b>　　%高調(diào)制頻率</b></p><p>　　y4=cos(n*0.8*pi);</p><p

37、>　　y5=y1.*y4';</p><p>　　f5=fft(y5);</p><p><b>　　%圖形顯示</b></p><p><b>　　figure</b></p><p>　　subplot(2,3,1),plot(t,y1)</p><p>　

38、　title('原信號時域圖')</p><p>　　xlabel('t')</p><p>　　ylabel('y1')</p><p>　　subplot(2,3,4),plot(w,abs(f1)*2/N)</p><p>　　title('原信號頻域圖')</p&g

39、t;<p>　　xlabel('w')</p><p>　　ylabel('Y1(w)')</p><p>　　subplot(2,3,2),plot(t,y3)</p><p>　　title('低調(diào)制后的時域圖')</p><p>　　xlabel('t')&l

40、t;/p><p>　　ylabel('y3')</p><p>　　subplot(2,3,5),plot(w,abs(f3)*2/N)</p><p>　　title('低調(diào)制后的頻域圖')</p><p>　　xlabel('w')</p><p>　　ylabel(&#

41、39;Y3(w)')</p><p>　　subplot(2,3,3),plot(t,y5)</p><p>　　title('高調(diào)制后的時域圖')</p><p>　　xlabel('t')</p><p>　　ylabel('y5')</p><p>　　su

42、bplot(2,3,6),plot(w,abs(f5)*2/N)</p><p>　　title('高調(diào)制后的頻域圖')</p><p>　　xlabel('w')</p><p>　　ylabel('Y5(w)')</p><p><b>　　%播放聲音</b></

43、p><p>　　wavplay(y1,fs)%播放調(diào)制前的原信號；</p><p>　　wavplay(y3,fs)%播放低調(diào)制后的信號；</p><p>　　wavplay(y5,fs)%播放高調(diào)制后的信號</p><p>　　【程序運行結果如下圖】：</p><p>　　分析：通過觀察原音樂信號的頻譜圖可知音樂信號的頻

44、率上限是0.5pi，但是為了方便以后的計算，在此將0.3pi后的信號舍去，即默認音樂信號的頻率上限是0.3pi。</p><p>　　取高頻調(diào)制頻率0.5pi，低頻調(diào)制頻率0.3pi對音樂信號進行調(diào)制。通過觀察不同調(diào)制頻率下的頻譜圖可以發(fā)現(xiàn)高頻調(diào)制的音樂信號頻譜發(fā)生了混疊，而采用低頻調(diào)制的音樂信號頻譜并未發(fā)生混疊。這是因為當采用高頻調(diào)制（0.5pi）時，頻譜被搬移到（2*n+1）*0,5pi，n=0.1.2.3…

45、..附近，此時高頻調(diào)制頻率高于原信號的頻率上限，故發(fā)生了頻譜混疊。同理，當采用低頻調(diào)制（0.3pi）時，未發(fā)生頻譜混疊。</p><p>　　播放不同調(diào)制頻率下的音樂信號，可以發(fā)現(xiàn)當采用低頻調(diào)制時，音樂信號比原信號的聲音低了很多，但是沒有雜音；采用高頻調(diào)制時，音樂信號比原信號的聲音低了很多的同時還伴隨有雜音。這是因為低頻調(diào)制沒有發(fā)生混疊，調(diào)制后的音樂信號頻譜幅度為原音樂信號的1/2，而高頻調(diào)制發(fā)生了混疊。<

46、/p><p>　　4、AM調(diào)制音樂信號的同步解調(diào)</p><p>　　① 設計巴特沃斯IIR濾波器完成同步解調(diào)；觀察濾波器頻率響應曲線；</p><p>　?、?用窗函數(shù)法設計FIR濾波器完成同步解調(diào)，觀察濾波器頻率響應曲線；（要求：分別使用矩形窗和布萊克曼窗，進行比較）； </p><p>　?、?輸出解調(diào)信號的波形和頻譜，觀察現(xiàn)象，給出理論解

47、釋； </p><p>　　④ 播放解調(diào)音樂信號，比較不同濾波器下的聲音，解釋現(xiàn)象。</p><p><b>　　理論基礎:</b></p><p>　　有一條信號f(t)恢復原始信號g(t)的過程稱為解調(diào)。這里，cos(w0*t)信號是接收端的本地載波信號，它與發(fā)送端的載波同頻同相。f(t)與cos(w0*t)相乘的結果使頻譜F(W)向左、右

48、分別移動+w0、-w0（并乘以系數(shù)1/2），得到 </p><p>　　g0(t)=1/2*g(t)+1/2*g(t)*cos(w0*t)</p><p>　　G0(w)=1/2*G(w)+1/4*[G(w-2*w0)+G(w+2*w0)]</p><p>　　再利用一個低通濾波器，濾除在頻率為2*w0附近的分量，即可取出g(t),完成解調(diào)。</p>

49、<p>　　【Matlab程序如下】：</p><p>　　%理想低通濾波器沖擊響應函數(shù)</p><p>　　function hd=ideal(N,wc)</p><p>　　for n=0:N-1</p><p>　　if n==(N-1)/2</p><p>　　hd(n+1)=wc/pi;<

50、/p><p>　　else hd(n+1)=sin(wc*(n-(N-1)/2))/(pi*(n-(N-1)/2));</p><p><b>　　end</b></p><p><b>　　end</b></p><p>　?。▽⑸鲜龀绦虮４鏋閕deal.m，但不能運行。）</p>&l

51、t;p><b>　　Clc</b></p><p><b>　　close all</b></p><p><b>　　clear all</b></p><p>　　[y,fs]=wavread('一生有你');</p><p>　　y1=y(:,1);

52、</p><p>　　N0=length(y);</p><p>　　n=0:(N0-1);</p><p>　　w=2*n/(N0-1);</p><p><b>　　t=n/fs;</b></p><p>　　f1=fft(y1);</p><p><b>　

53、　%調(diào)制頻率濾波</b></p><p>　　y2=cos(n*0.4*pi);%調(diào)制頻率為0.4*pi;</p><p>　　y3=y1.*y2';%調(diào)制后的信號;</p><p>　　y4=y3.*y2';%調(diào)制解調(diào)后的信號;</p><p>　　%對調(diào)制頻率巴特沃斯IIR濾波器的設計</p>&

54、lt;p><b>　　Wp=0.2;</b></p><p><b>　　Ws=0.3;</b></p><p><b>　　Rp=1;</b></p><p><b>　　Rs=20;</b></p><p>　　W=0:0.001*pi:1*pi

55、;</p><p>　　[N,Wc]=buttord(Wp,Ws,Rp,Rs);</p><p>　　[b,a]=butter(N,Wc);</p><p>　　[H,W]=freqz(b,a,W);</p><p>　　plot(W/pi,abs(H))</p><p>　　title('巴特沃斯濾波器的頻域

56、圖')</p><p>　　xlabel('w')</p><p>　　ylabel('H(w)')</p><p>　　y5=filter(b,a,y4)*2;%對調(diào)制解調(diào)后的信號進行濾波;</p><p>　　f3=fft(y3);%對調(diào)制信號做fft變換;</p><p>

57、　　f4=fft(y4);%對調(diào)制解調(diào)后的信號做fft變換;</p><p>　　f5=fft(y5);%對調(diào)制解調(diào)濾波后的信號做fft變換;</p><p><b>　　figure</b></p><p>　　subplot(4,1,1),plot(w,abs(f1)*2/N0)%調(diào)制前的原信號頻域圖；</p><p&g

58、t;　　title('原信號頻域圖')</p><p>　　xlabel('w')</p><p>　　ylabel('Y1(w)')</p><p>　　subplot(4,1,2),plot(w,abs(f3)*2/N0)%高調(diào)制后的信號頻域圖；</p><p>　　title('調(diào)

59、制后的頻域圖')</p><p>　　xlabel('w')</p><p>　　ylabel('Y3(w)')</p><p>　　subplot(4,1,3),plot(w,abs(f4)*2/N0)</p><p>　　title('調(diào)制解調(diào)后的頻域圖')</p>&

60、lt;p>　　xlabel('w')</p><p>　　ylabel('Y4(w)')</p><p>　　subplot(4,1,4),plot(w,abs(f5)*2/N0) %解調(diào)濾波后的信號頻域圖；</p><p>　　title('解調(diào)濾波后的頻域圖')</p><p>　　

61、xlabel('w')</p><p>　　ylabel('Y5(w)')</p><p>　　%用窗函數(shù)設計FIR濾波器</p><p><b>　　N=33; </b></p><p>　　Wc1=0.4*pi;</p><p>　　hd=ideal(N,W

62、c1);</p><p>　　w1=boxcar(N);%矩形窗</p><p>　　w2=blackman(N);%布萊克曼窗</p><p>　　h1=hd.*w1';</p><p>　　h2=hd.*w2';</p><p><b>　　%求濾波器頻率響應</b></

63、p><p><b>　　M=1024;</b></p><p>　　W=2/M*[0:M-1];</p><p>　　fh1=fft(h1,M);</p><p>　　db1=-20*log10(abs(fh1(1)./(abs(fh1)+eps))); %幅度的分貝表示</p><p>　

64、　fh2=fft(h2,M);</p><p>　　db2=-20*log10(abs(fh2(1)./(abs(fh2)+eps))); %幅度的分貝表示</p><p><b>　　figure </b></p><p>　　subplot(2,2,1);plot(W,abs(fh1))</p><p>　

65、　title('矩形窗幅度')</p><p>　　subplot(2,2,2);plot(W,db1)</p><p>　　title('矩形窗分貝') </p><p>　　subplot(2,2,3);plot(W,abs(fh2))</p><p>　　title('布萊克曼窗幅度')&

66、lt;/p><p>　　subplot(2,2,4);plot(W,db2)</p><p>　　title('布萊克曼窗分貝')</p><p>　　y6=conv(y4,h1);%作卷積</p><p>　　f6=fft(y6);</p><p>　　y7=conv(y4,h2);</p>

67、<p>　　f7=fft(y7);</p><p>　　M1=length(y4)+length(h1)-1;</p><p>　　W1=2/M1*[0:M1-1];</p><p>　　M2=length(y4)+length(h2)-1;</p><p>　　W2=2/M2*[0:M2-1];</p><p

68、><b>　　figure</b></p><p>　　subplot(3,1,1),plot(w,abs(f1)*2/(length(y)))</p><p>　　title('原信號頻域圖')</p><p>　　xlabel('w')</p><p>　　ylabel('

69、;Y1(w)')</p><p>　　subplot(3,1,2),plot(W1,abs(f6)*2/M1)</p><p>　　title('用矩形窗濾波后的頻域圖')</p><p>　　xlabel('w')</p><p>　　ylabel('Y6(w)')</p>

70、<p>　　subplot(3,1,3),plot(W2,abs(f7)*2/M2)</p><p>　　title('用布萊克曼窗濾波后的頻域圖')</p><p>　　xlabel('w')</p><p>　　ylabel('Y7(w)')</p><p>　　wavplay

71、(y1,fs)</p><p>　　wavplay(y6,fs)</p><p>　　wavplay(y7,fs)</p><p><b>　　程序運行如下：</b></p><p>　　分析1：通過對比解調(diào)前后的音樂信號時域圖和頻譜圖可知，解調(diào)后的音樂信號時域和頻譜幅度變?yōu)樵瓉淼?/2，若想濾波后可以完全恢復原始音樂信

72、號，可以在傳輸函數(shù)添加系數(shù)2即可。</p><p>　　分析2：對FIR窗函數(shù)濾波器的設計參數(shù)進行選定：</p><p>　　矩形窗（boxcar）：濾波器的截止頻率wc=0.5pi，濾波器的階數(shù)N=33，阻帶衰減為25db；布萊克曼窗（blackman）：濾波器的截止頻率wc=0.5pi，濾波器的階數(shù)N=33，阻帶衰減為80db。</p><p>　　比較采用矩形

73、窗和布萊克曼窗的頻率特性圖可以看出：最小阻帶衰減只由窗形狀決定，而不受階數(shù)N的影響；而過渡帶的寬度窗的形狀有關；同時，布萊克曼窗的過濾帶寬、旁瓣峰值和主瓣寬度均大于矩形窗的過濾帶寬、旁瓣峰值和主瓣寬度。</p><p>　　分析3：下圖為經(jīng)過加矩形窗和布萊克曼窗的濾波器濾波后的音樂信號的時域圖和頻譜圖。比較兩種窗函數(shù)濾波下的頻譜圖可知加矩形窗的濾波器的濾波效果明顯沒有加布萊克曼窗的濾波器濾波效果好，原因如分析2所

74、示：布萊克曼窗的最小阻帶衰減在同樣階數(shù)的前提條件下要比矩形窗的最小阻帶衰減大，所以濾波更徹底，恢復的原信號也就會比矩形窗好。若想提高矩形窗的濾波效果可采取增加階數(shù)的方法，但這種情況下會增加不必要的開支，所以在選定濾波器的時候要權衡各方面的條件，力爭達到一個各方都滿意的結果。</p><p>　　分析4：播放IIR濾波器和FIR濾波器濾波后的音樂信號，我們可以聽到音樂信號的音調(diào)變低了，這是因為高頻部分被濾掉了。同時

75、還可以聽到IIR濾波器濾波后的音樂信號沒有FIR濾波器濾波后的音樂信號效果好，這是因為FIR濾波器是線性相位的，而IIR濾波器是非線性相位；矩形窗沒有布萊克曼窗的濾波效果好。</p><p>　　5、音樂信號的濾波去噪</p><p>　?、俳o原始音樂信號疊加幅度為0.05，頻率為3kHz、 5kHz、8kHz的三余弦混合噪聲，觀察噪聲頻譜以及加噪后音樂信號的音譜和頻譜，并播放音樂，感受噪

76、聲對音樂信號的影響；</p><p>　?、诮o原始音樂信號疊加幅度為0.5的隨機白噪聲（可用rand語句產(chǎn)生），觀察噪聲頻譜以及加噪后音樂信號的音譜和頻譜，并播放音樂，感受噪聲對音樂信號的影響；</p><p>　?、鄹鶕?jù)步驟①、②觀察到的頻譜，選擇合適指標設計濾波器進行濾波去噪，觀察去噪后信號音譜和頻譜，并播放音樂，解釋現(xiàn)象。 </p><p><b>

77、　　理論基礎：</b></p><p>　　Rand函數(shù)介紹：rand函數(shù)產(chǎn)生由在（0，1）之間均勻分布的隨機數(shù)組組成的 數(shù)組。Y = rand(n) 返回一個n x n的隨機矩陣。如果n不是數(shù)量，則返回錯誤信息。</p><p>　　【Matlab程序如下】：</p><p><b>　　Clc</b&g

78、t;</p><p><b>　　Close all</b></p><p><b>　　clear all</b></p><p>　　[y,fs]=wavread('一生有你');</p><p>　　y1=y(:,1);</p><p>　　f1=fft

79、(y1);</p><p>　　N0=length(y);</p><p>　　n=0:(N0-1);</p><p><b>　　w=2*n/N0;</b></p><p><b>　　t=n/fs;</b></p><p>　　y2=0.05*cos(2*pi*3000*

80、t)+0.05*cos(2*pi*5000*t)+0.05*cos(2*pi*8000*t);%三余弦函數(shù)</p><p>　　f2=fft(y2);</p><p>　　y3=y1+y2';</p><p>　　f3=fft(y3);</p><p><b>　　figure</b></p>&l

81、t;p>　　subplot(3,2,1),plot(t(1:1000),y1(1:1000))</p><p>　　title('原信號時域圖')</p><p>　　xlabel('t')</p><p>　　ylabel('y1')</p><p>　　subplot(3,2,2),p

82、lot(w,abs(f1)*2/N0)</p><p>　　title('原信號頻譜')</p><p>　　xlabel('w')</p><p>　　ylabel('Y1(w)')</p><p>　　subplot(3,2,3),plot(t(1:200),y2(1:200))</p

83、><p>　　title('三余弦信號時域圖')</p><p>　　xlabel('t')</p><p>　　ylabel('y2')</p><p>　　subplot(3,2,4),plot(w,abs(f2)*2/N0)</p><p>　　title('三

84、余弦信號頻譜')</p><p>　　xlabel('w')</p><p>　　ylabel('Y2(w)')</p><p>　　subplot(3,2,5),plot(t(1:1000),y3(1:1000))</p><p>　　title('加噪后信號的時域')</p&g

85、t;<p>　　xlabel('t')</p><p>　　ylabel('y2')</p><p>　　subplot(3,2,6),plot(w,abs(f3)*2/N0)</p><p>　　title('加噪后信號的的頻譜')</p><p>　　xlabel('w

86、')</p><p>　　ylabel('Y3(w)')</p><p>　　y4=rand(N0,1)-0.5;%白噪聲的產(chǎn)生算法</p><p><b>　　y5=y1+y4;</b></p><p>　　f4=fft(y4);</p><p>　　f5=fft(y5)

87、;</p><p><b>　　figure</b></p><p>　　subplot(2,2,1),plot(t(1:1000),y4(1:1000))</p><p>　　title('隨機白噪聲時域圖')</p><p>　　xlabel('t')</p><p

88、>　　ylabel('y4')</p><p>　　subplot(2,2,3),plot(t(1:1000),y5(1:1000))</p><p>　　title('帶有隨機白噪聲的信號時域')</p><p>　　xlabel('t')</p><p>　　ylabel('y

89、5')</p><p>　　subplot(2,2,2),plot(w,abs(f4)*2/N0)</p><p>　　title('隨機白噪聲頻譜')</p><p>　　xlabel('w')</p><p>　　ylabel('Y4(w)')</p><p>

90、;　　subplot(2,2,4),plot(w,abs(f5)*2/N0)</p><p>　　title('帶有隨機白噪聲的信號頻譜')</p><p>　　xlabel('w')</p><p>　　ylabel('Y5(w)')</p><p>　　%巴特沃斯濾波器的設計</p&g

91、t;<p><b>　　Wp=0.08;</b></p><p>　　Ws=0.12;%不宜過狹窄 不然會出錯</p><p><b>　　Rp=1;</b></p><p>　　Rs=30 ; %指標輸入</p><p>　　[N,Wc]=buttord(Wp,Ws,Rp,Rs);

92、</p><p>　　[b,a]=butter(N,Wc);</p><p>　　[H,W]=freqz(b,a);</p><p><b>　　figure</b></p><p>　　plot(W/pi,abs(H))</p><p>　　title('低巴特沃斯濾波器的頻域圖'

93、;)</p><p>　　xlabel('w')</p><p>　　ylabel('H(w)')</p><p>　　y6=filter(b,a,y3);</p><p>　　y7=filter(b,a,y5);</p><p>　　f6=fft(y6);</p><

94、;p>　　f7=fft(y7);</p><p><b>　　figure</b></p><p>　　subplot(3,1,1),plot(w,abs(f1)*2/N0)</p><p>　　title('原信號頻譜')</p><p>　　xlabel('w')</p&g

95、t;<p>　　ylabel('Y1(w)')</p><p>　　subplot(3,1,2),plot(w,abs(f6)*2/N0)</p><p>　　title('濾除三余弦信號后的頻譜')</p><p>　　xlabel('w')</p><p>　　ylabel(&

96、#39;Y6(w)')</p><p>　　subplot(3,1,3),plot(w,abs(f7)*2/N0)</p><p>　　title('濾除白噪聲信號后的頻譜')</p><p>　　xlabel('w')</p><p>　　ylabel('Y7(w)')</p&g

97、t;<p>　　%wavplay(y3,fs) 加三余弦后 一直有個比較尖銳刺耳的聲音存在 像蜂鳴聲</p><p>　　%wavplay(y6,fs) 濾除后跟原音樂旋律一樣 只是聲調(diào)低沉些</p><p>　　%wavplay(y5,fs) 加白噪聲后 就像電視機有雪花收不到臺時的那種嘈噪聲 </p><p>　　%wavplay(y7,fs

98、) 濾除后噪聲仍存在但聲音低了些</p><p>　　【程序運行結果如下圖】：</p><p>　　分析1：上圖為三余弦信號的時域圖和頻譜圖以及加三余弦噪聲后的音樂信號</p><p>　　時域圖和頻譜圖。通過計算可知：三余弦混合噪聲信號的三條頻譜線分別加在頻率為3kz，5kz和8kz處。而且通過播放加三余弦噪聲的音樂信號可知，同原音樂信號相比增加了尖銳的刺耳聲，

99、這同預想的信號一致。</p><p>　　分析2：下圖為隨機白噪聲信號的時域圖和頻譜圖以及加隨機白噪聲信號后的</p><p>　　音樂信號時域圖和頻譜圖。而且通過播放加隨即白噪聲信號的音樂信號可知，同原音樂信號相比增加了嗞嗞的電波聲，這同預想的信號一致。</p><p>　　分析3：隨機白噪聲去噪：由于隨即白噪聲信號是加在了整個時域和頻域內(nèi)，為了濾波的更徹底，所以

100、選取較小的通帶截止頻率和阻帶截止頻率，并選取較大的組帶最小衰減。wp=0.08pi為通帶截止頻率，ws=0.12pi為其阻帶截止頻率，rp=1db為其通帶最大衰減,，rs=30db為阻帶最小衰減。</p><p>　　分析4：觀察下圖三余弦噪聲去噪后的時域圖和頻譜圖，以及隨機白噪聲去噪后的時域圖和頻譜圖可知：三余弦噪聲被完全濾掉；隨機白噪聲由于是加在原因樂信號的整個時域和頻域內(nèi)，所以濾波后的音樂信號仍然殘留有少量

101、的噪聲信號，這是避免不了的。把去噪后的音樂信號時域圖和頻譜圖和原始信號相比發(fā)現(xiàn)并無大的差別，所以在此可以暫定為濾波成功。</p><p>　　音樂信號的幅頻特性和相頻分析</p><p>　　設計低通濾波器，濾除原始信號高頻信息；設計高通濾波器，濾除原始信號低頻信息；選取兩端不同的音樂信號，分別將其幅度譜與相位譜交叉組合構成新的音樂信號，播放并比較組合后的音樂與原始信號，感受相頻信息對音樂

102、信號的影響</p><p>　　【matlab的程序】:</p><p><b>　　clc</b></p><p><b>　　close all</b></p><p><b>　　clear all</b></p><p>　　[y,fs]=wav

103、read('一生有你');</p><p>　　y1=y(:,1);</p><p>　　f1=fft(y1);</p><p>　　N0=length(y);</p><p>　　n=0:(N0-1);</p><p><b>　　w=2*n/N0;</b></p>

104、<p><b>　　t=n/fs;</b></p><p>　　%低通巴特沃斯濾波器的設計</p><p><b>　　Wp1=0.1;</b></p><p><b>　　Ws1=0.2;</b></p><p><b>　　Rp=1;</b>

105、;</p><p><b>　　Rs=20;</b></p><p>　　W=0:0.001*pi:0.5*pi;</p><p>　　[N1,Wc1]=buttord(Wp1,Ws1,Rp,Rs);</p><p>　　[b1,a1]=butter(N1,Wc1);</p><p>　　[H1,

106、W]=freqz(b1,a1,W);</p><p>　　y2=filter(b1,a1,y1);</p><p>　　f2=fft(y2);</p><p>　　%高通巴特沃斯濾波器的設計</p><p><b>　　Wc2=0.3;</b></p><p><b>　　Ws2=0.2

107、;</b></p><p><b>　　Rc=1;</b></p><p><b>　　Rs=20;</b></p><p>　　W=0:0.001*pi:0.5*pi;</p><p>　　[N2,Wc2]=buttord(Wc2,Ws2,Rc,Rs);</p><

108、p>　　[b2,a2]=butter(N2,Wc2,'high');</p><p>　　[H2,W]=freqz(b2,a2,W);</p><p>　　y3=filter(b2,a2,y1);</p><p>　　f3=fft(y3);</p><p><b>　　figure</b></

109、p><p>　　subplot(4,2,1),plot(t(1:1000),y1(1:1000))</p><p>　　title('原信號時域圖')</p><p>　　xlabel('t')</p><p>　　ylabel('y1')</p><p>　　subplot

110、(4,2,2),plot(w,abs(f1)*2/N0)</p><p>　　title('原信號頻譜')</p><p>　　xlabel('w')</p><p>　　ylabel('Y1(w)')</p><p>　　subplot(4,2,3),plot(W/pi,abs(H1))<

111、;/p><p>　　title('低通巴特沃斯濾波器的頻域圖')</p><p>　　xlabel('w')</p><p>　　ylabel('H(w)')</p><p>　　subplot(4,2,5),plot(t(1:1000),y2(1:1000))</p><p&g

112、t;　　title('濾除高頻后的時域圖')</p><p>　　xlabel('t')</p><p>　　ylabel('y2')</p><p>　　subplot(4,2,7),plot(w,abs(f2)*2/N0)</p><p>　　title('濾除高頻后的頻譜'

113、)</p><p>　　xlabel('w')</p><p>　　ylabel('Y2(w)')</p><p>　　subplot(4,2,4),plot(W/pi,abs(H2))</p><p>　　title('高通巴特沃斯濾波器的頻域圖')</p><p>　

114、　xlabel('w')</p><p>　　ylabel('H2(w)')</p><p>　　subplot(4,2,6),plot(t(1:1000),y3(1:1000))</p><p>　　title('濾除低頻后的時域圖')</p><p>　　xlabel('t'

115、)</p><p>　　ylabel('y3')</p><p>　　subplot(4,2,8),plot(w,abs(f3)*2/N0)</p><p>　　title('濾除低頻后的頻譜')</p><p>　　xlabel('w')</p><p>　　ylabe

116、l('Y3(w)')</p><p>　　%wavplay(y1,fs) %原信號聲音</p><p>　　%wavplay(y2,fs) %濾除高頻后的信號聲音 聲音有點低沉</p><p>　　%wavplay(y3,fs) %濾除低頻后信號的聲音 聲音有點尖銳</p><p><b>　　%交叉組合&

117、lt;/b></p><p>　　[x1,fs]=wavread('一生有你');</p><p>　　[x2,fs]=wavread('天路');</p><p>　　y1=x1(:,1);</p><p>　　y2=x2(:,1);</p><p>　　y3=y2(1:leng

118、th(y1));%截取長度使得長度相等</p><p>　　N=length(y1);</p><p>　　f1=fft(y1);</p><p>　　f2=fft(y3);</p><p>　　Fy1=abs(f1);%幅度</p><p>　　Ay1=angle(f1);</p><p>　

119、　Fy3=abs(f2);</p><p>　　Ay3=angle(f2);%角度</p><p>　　F1=Fy1.*exp(j.*Ay3);</p><p>　　F2=Fy3.*exp(j.*Ay1);</p><p>　　y11=ifft(F1);</p><p>　　y22=ifft(F2);</p>

120、;<p>　　w=2/N*[0:N-1];</p><p>　　wavplay(x1,fs);</p><p>　　wavplay(real(y11),fs);%聲音大致能聽到是第二個的聲音 但背影很嘈雜</p><p>　　wavplay(x2,fs);</p><p>　　wavplay(real(y22),fs);%聲音

121、大致能聽到是第一個的聲音 但背影很嘈雜</p><p>　　%結論 相位角決定不同的聲音 可用來區(qū)分聲音 音調(diào)即聲音大小取決于幅度且成正相關</p><p>　　Matlab的運行結果：</p><p><b>　　四、問題討論</b></p><p>　　1、IIR數(shù)字濾波器和FIR數(shù)字濾波器的比較</p&

122、gt;<p>　　答：IIR濾波器：相位一般是非線性的，不一定穩(wěn)定，不能使用fft做快速卷積，一定是遞歸結構，對頻率分量的選擇性好（零極點可同時起作用），相同性能下階次較低，有噪聲反饋，噪聲大，運算誤差大，有可能出現(xiàn)極限環(huán)振蕩，設計時有大量圖表可查，方便簡單，主要用于設計分段常數(shù)的標準低通，帶通，高通，帶阻和全通濾波器。</p><p>　　FIR濾波器：相位可以做到嚴格線性，一定是穩(wěn)定的，信號通過

123、系統(tǒng)可采用快速卷積，主要是非遞歸結構，也可含遞歸結構，選擇性差，相同性能下階次高，噪聲小，運算誤差小，不會出現(xiàn)極限環(huán)振蕩，可設計正交變換器、微分器、線性預測器、線性調(diào)頻器等各種網(wǎng)絡，適用范圍廣。</p><p>　　其中，F(xiàn)IR濾波器的最大好處是穩(wěn)定、線性相位和廣泛的適用范圍，而它的最大</p><p>　　缺點是階數(shù)高，從而帶來時延大、存儲單元多等問題。例如用頻率抽樣法設計阻帶衰耗為－2

124、0dB的FIR DF需33階，用雙線性法設計同樣指標的切比雪夫IIR DF僅需4～5階。因此，在一些對時延有嚴格限制的場合就不得不考慮用IIR濾波器。語音信號對相位的非線性不很敏感。數(shù)據(jù)和圖象信號則往往對濾波器提出線性相位的要求，這就是為什么FIR用得越來越廣的原因?？傊?，IIR和FIR各有特點，在應用時要根據(jù)各方面的指標，綜合考慮加以選擇。</p><p>　　2、音樂信號的音調(diào)與信號的什么特征有關？</

125、p><p>　　答：音調(diào)主要由聲音的頻率決定。對一定強度的純音，音調(diào)隨頻率的升降而升降；對一定頻率的純音、低頻純音的音調(diào)隨響度增加而下降，高頻純音的音調(diào)卻隨響度增加而上升。音調(diào)的高低還與發(fā)聲體的結構有關，因為發(fā)聲體的結構影響了聲音的頻率。所以音樂信號的音調(diào)與信號的頻率有關。音調(diào)還與聲音持續(xù)的時間長短有關。非常短促(毫秒量級或更短)的純音，只能聽到像打擊或彈指那樣的“喀嚓”一響，感覺不出音調(diào)。持續(xù)時間從10 毫秒增加到

126、50 毫秒，聽起來覺得音調(diào)是由低到高連續(xù)變化的。超過50 毫秒，音調(diào)就穩(wěn)定不變了。</p><p>　　3、音樂信號的音色與信號的什么特征有關？</p><p>　　答：音樂信號的音色是指對聲音音質(zhì)的感覺特性，也是一種聲音區(qū)別于另一種聲音的特征品質(zhì)。音色的不同取決于不同的泛音，每一種樂器、不同的人以及所有能發(fā)聲的物體發(fā)出的聲音，除了一個基音還有許多不同頻率的泛音伴隨，正是這些泛音決定了其不

127、同的音色，使人能辨別出是不同的樂器甚至不同的人發(fā)出的聲音。而泛音的不同指的是頻譜的不同，所以音樂信號的音色與信號的頻譜或頻率特性有關。</p><p>　　4、兩種不同音色的音樂信號疊加混疊后，為何人耳還可以分辨？</p><p>　　答：人耳聽覺系統(tǒng)實際上就是一個音頻信號處理器，可以完成對聲音信號的傳輸、轉換以及綜合處理的能力。而音樂信號的音色與音樂信號的頻譜分布有關，不同音色的音樂信號

128、的頻譜當然不會相同，當混疊后的音樂信號進入人耳的聽覺系統(tǒng)后，人耳通過對混疊信號的頻譜處理，可以將具有不同頻譜的聲音信號分離開，也就達到了分辨不同音色的音樂信號的目的。</p><p>　　5、音樂信號的幅度與相位特征對信號有哪些影響？</p><p>　　答：音樂信號的相位發(fā)生改變會引起信號時域特性發(fā)生改變；音樂信號的幅度發(fā)生變化會影響信號的頻譜幅度。</p><p&g

129、t;<b>　　五、心得體會</b></p><p>　　這次實驗是以上學期所學的數(shù)字信號理論為依據(jù)，用matlab程序?qū)崿F(xiàn)，通過多次上機編程、運行、調(diào)試，學到了很多東西。首先，為了完成這次的課程設計，需要從新將上學期遺忘的知識給撿起來，而且還需要去圖書館查閱有關matlab編程的相關知識，學習有關函數(shù)的用法；其次，在做實驗的過程中會遇到很多的困難，在困難面前不要放棄，只要有信心和耐心，堅持