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

1、長(zhǎng)春理工大學(xué)電工電子實(shí)驗(yàn)教學(xué)中心學(xué)生實(shí)驗(yàn)報(bào)告2014 2015 學(xué)年第一 學(xué)期實(shí)驗(yàn)課程數(shù)字信號(hào)處理實(shí)驗(yàn)實(shí)驗(yàn)地點(diǎn) 東4教學(xué)樓444實(shí)驗(yàn)室學(xué)院電子信息工稈學(xué)院專業(yè)通信工稈學(xué)號(hào)120421101姓緒楊杰實(shí)驗(yàn)項(xiàng)目音頻信號(hào)釆集、分析及處理實(shí)驗(yàn)時(shí)間第十四周星期四9、10節(jié)實(shí)驗(yàn)臺(tái)號(hào)操作成績(jī)報(bào)告成績(jī)?nèi)龑?shí)驗(yàn)?zāi)康?. 以音頻信號(hào)為例，熟悉模擬信號(hào)數(shù)字處理過程，進(jìn)一步理解數(shù)字信號(hào)處理概念。2. 掌握運(yùn)用matlab實(shí)現(xiàn)對(duì)音頻信號(hào)的時(shí)頻分析方法；3. 初步掌握數(shù)字音頻信號(hào)合成的方法。4. 掌握運(yùn)用matlab設(shè)計(jì)iir和fir濾波系統(tǒng)的方法；5. 掌握運(yùn)用matlab實(shí)現(xiàn)對(duì)加噪音頻信號(hào)進(jìn)行去噪濾波的方法。鍛煉學(xué)生運(yùn)用

2、所學(xué)知識(shí)獨(dú)立分析問題解決問題的能力，培養(yǎng)學(xué)生創(chuàng)新能力。二、實(shí)驗(yàn)儀器計(jì)算機(jī)、matlab軟件三、實(shí)驗(yàn)性質(zhì)綜合分析、設(shè)計(jì)性實(shí)驗(yàn)四、實(shí)驗(yàn)任務(wù)實(shí)驗(yàn)內(nèi)容一：windows系統(tǒng)中的“ding”音頻信號(hào)的采集、分析、合成1. 音頻信號(hào)的采集編寫matlab程序，采集windows系統(tǒng)中的“ ding”聲，得到*.wav音頻文件，而后實(shí)現(xiàn)音頻信 號(hào)回放。2. 音頻信號(hào)的頻譜分析運(yùn)用matlab軟件實(shí)現(xiàn)對(duì)音頻信號(hào)的時(shí)域分析和頻域分析，并打印相應(yīng)的圖形，完成在實(shí)驗(yàn)報(bào) 告中。注意：此音頻信號(hào)的頻譜包含兩條主要譜線，在進(jìn)行頻譜分析時(shí)，注意頻譜的完整性，利用matlab實(shí)現(xiàn)對(duì)兩條主要譜線的定位并計(jì)算譜線所對(duì)應(yīng)的模擬頻率

3、。3. 音頻信號(hào)的分解和合成運(yùn)用matlab軟件實(shí)現(xiàn)音頻信號(hào)的分解與合成，將音頻信號(hào)的頻譜中兩部分頻譜成分進(jìn)行分解，分別繪制出分解后的兩個(gè)信號(hào)的頻譜圖；然后將分解后的兩個(gè)信號(hào)再合成為一個(gè)新的信號(hào)，將合成后的新信號(hào)的時(shí)域、頻域圖與原來的信號(hào)時(shí)域、頻域圖相比較，繪制岀對(duì)比效果圖。4. 音頻信號(hào)的回放運(yùn)用matlab軟件實(shí)現(xiàn)音頻信號(hào)的回放，將合成后的新信號(hào)和原音頻信號(hào)分別進(jìn)行回放，對(duì)比兩個(gè)信號(hào)的聲音效果。5. 音頻信號(hào)分段傅里葉分析（選作）分析對(duì)一般音頻.wav信號(hào)進(jìn)行一次性傅里葉分析時(shí)存在的主要問題，利用分段傅里葉變換對(duì)該音頻信號(hào)重新分析并合成。對(duì)比一次傅里葉分析結(jié)果并進(jìn)行總結(jié)。實(shí)驗(yàn)內(nèi)容二：任意音

4、頻信號(hào)的時(shí)域和頻域分析及數(shù)字濾波器設(shè)計(jì)1. 音頻信號(hào)的采集音頻信號(hào)的采集可以通過windows自帶的錄音機(jī)也可以用專用的錄制軟件錄制一段音頻信號(hào)（盡量保證無噪音、干擾?。部梢灾苯訌?fù)制一段音頻信號(hào)（時(shí)間為1s）,但必須保證音頻信號(hào)保存為.wav的文件。2. 音頻信號(hào)的時(shí)域、頻域分析運(yùn)用matlab軟件實(shí)現(xiàn)對(duì)音頻信號(hào)的打開操作以及時(shí)域分析和頻域分析，并畫出相應(yīng)的圖形，打印在實(shí)驗(yàn)報(bào)告中。3. 引入干擾信號(hào)在原有的音頻信號(hào)上，疊加一個(gè)頻率為100khz的正弦波干擾信號(hào)（幅度自定，可根據(jù)音頻信號(hào)幅度情況而定）。4. 數(shù)字濾波器設(shè)計(jì)運(yùn)用matlab實(shí)現(xiàn)iir或fir濾波系統(tǒng)，要求加入干擾的音頻信號(hào)經(jīng)過

5、濾波系統(tǒng)后，能夠?yàn)V除 100khz的干擾信號(hào)，同時(shí)保留原有的音頻信號(hào)，要求繪制出濾波系統(tǒng)的頻譜圖，同時(shí)繪制出濾波 前后音頻信號(hào)的時(shí)域、頻域?qū)Ρ葓D。實(shí)驗(yàn)內(nèi)容三：編寫 fft算法程序（此部分為選作內(nèi)容）分析fft算法實(shí)現(xiàn)過程，運(yùn)用 matlab軟件實(shí)現(xiàn)fft函數(shù)功能，要求運(yùn)用自行編寫的matlab程序?qū)崿F(xiàn)對(duì)固定信號(hào)（自定義信號(hào)）的頻譜分析，并繪制出頻譜圖；然后將自己設(shè)計(jì)的fft算法程序與fft函數(shù)實(shí)現(xiàn)的頻譜圖相比較，分析兩個(gè)程序結(jié)果的差別及原因。五、實(shí)驗(yàn)任務(wù)分解及思路實(shí)驗(yàn)內(nèi)容一：windows系統(tǒng)中的“ding”音頻信號(hào)的釆集、分析、合成通過對(duì)任務(wù)內(nèi)容要求的分析，建立程序框圖如圖一（任務(wù)一程序框圖

6、）所示，詳細(xì)程序及結(jié)果見附錄。編程思想及思路如下：1、首先用audioread讀取音頻信號(hào)得到 signal序列及采樣頻率fs,由于音頻信號(hào)是雙聲道，這里只取第一個(gè)聲道進(jìn)行處理，由于原信號(hào)本身就較短，這里對(duì)長(zhǎng)度就不作處理。接著按實(shí)驗(yàn)要求進(jìn)行音頻回放，并繪制出音頻信號(hào)的時(shí)域波形。2、調(diào)用fft函數(shù)，對(duì)signal進(jìn)行快速傅里葉變換，并繪制出頻譜圖。在繪制頻譜圖時(shí)考慮到快 速傅里葉變換的對(duì)稱性，這里只取序列的前半部分進(jìn)行觀察分析。3、為了實(shí)現(xiàn)音頻信號(hào)的分解及合成，先對(duì)原信號(hào)的頻譜圖進(jìn)行觀察分析，發(fā)現(xiàn)原信號(hào)的主要能量集中在兩個(gè)頻率上，于是考慮用這兩個(gè)頻率的正弦信號(hào)合成該信號(hào)。這樣首先要找到 這兩個(gè)頻

7、率。為了求得這兩個(gè)頻率，我先調(diào)用findpeaks函數(shù)找到頻譜圖上的各個(gè)峰值點(diǎn)signal_pks及其對(duì)應(yīng)的位置iocs ,然后用sort對(duì)峰值點(diǎn)進(jìn)行排序，找到最大的兩個(gè)值， 接著用find函數(shù)找到這兩個(gè)峰值點(diǎn)在iocs中的位置，這樣就從iocs中知道了最大兩個(gè)峰值點(diǎn)在原音頻信號(hào)中的位置，也就知道了對(duì)應(yīng)的頻率。接著生成這兩個(gè)頻率的正弦信號(hào)，幅度的比例等于頻譜圖上對(duì)應(yīng)頻率的幅值比例。然后將兩個(gè)正弦信號(hào)疊加，得到合成后的信號(hào)，繪制出時(shí)域波形，與原信號(hào)波形進(jìn)行比較，接著對(duì)兩個(gè)正弦信號(hào)進(jìn)行fft ,繪制岀他們的頻譜，然后對(duì)合成的信號(hào)進(jìn)行fft ,做出頻譜圖與原信號(hào)的頻譜圖進(jìn)行比較。4、用sound函數(shù)

8、進(jìn)行原信號(hào)和合成信號(hào)的回放，比較差異。5、要對(duì)音頻信號(hào)進(jìn)行分段傅里葉分析，首先就要確定每段的長(zhǎng)度?？紤]計(jì)算機(jī)進(jìn)行傅里葉變換時(shí)點(diǎn)數(shù)為時(shí)，為便于計(jì)算，這里就將每段長(zhǎng)度暫定為1024,在for循環(huán)中對(duì)每段信號(hào)用fft函數(shù)進(jìn)行快速傅里葉變換，然后再用ifft對(duì)每段快速傅里葉變換后的序列進(jìn)行快速傅里葉逆變換并合成，然后進(jìn)行音頻回放，與原信號(hào)進(jìn)行比較。尋找頻譜最高護(hù)個(gè)峰值點(diǎn)找到峰值點(diǎn)苓應(yīng)的頻率等于原信號(hào)最大幅值做合成信號(hào)時(shí)域波 些原罕出交 做合成信號(hào)丁域旅丈圖 士故兩旺弦罕的頻譜做合成信號(hào)的頻譜與原信嗎華i普比孫原信號(hào)及合成信號(hào)的回護(hù)較原信號(hào)的分段傅里圖一任務(wù)一程序框圖實(shí)驗(yàn)內(nèi)容二：任意音頻信號(hào)的時(shí)域和頻域

9、分析及數(shù)字濾波器設(shè)計(jì)通過對(duì)任務(wù)內(nèi)容要求的分析，建立程序框圖如圖二（任務(wù)二程序框圖）所示，詳細(xì)程序及結(jié)果見附錄。編程思想及思路如下：1、首先用audioread讀取音頻信號(hào)signal及采樣頻率fs,同實(shí)驗(yàn)內(nèi)容樣只取第一聲道信號(hào)2的數(shù)據(jù)進(jìn)行處理。進(jìn)行分析處理，由于此處選擇的音頻文件較長(zhǎng)，這里截取前2、同實(shí)驗(yàn)任務(wù)一一樣，做出信號(hào)的時(shí)域波形，及 fft變換后的頻譜圖。圖二任務(wù)二程序框圖3、這里要求引入100khz的正弦干擾信號(hào)，由于由 1中已得到fs為22050h乙 根據(jù)奈奎斯特頻率釆樣定理，采樣頻率必須大于等于原信號(hào)最高頻率的2倍，所以必須將原信號(hào)的采樣頻率提高到200khz以上才能引入100kh

10、z的噪聲，所以這里考慮用一階線型插值interpl將原信號(hào)的采樣頻率提高到220500h乙這樣就可以引入 100khz噪聲。做出提高采樣頻率后的信號(hào)的時(shí)域波形和頻譜圖，確認(rèn)信號(hào)并沒有發(fā)生變化。接著生成100khz的正弦信號(hào)，根據(jù)2中做出的信號(hào)的時(shí)域波形的幅度，這里取噪聲的幅值為0.5o將提高采樣頻率后的4、信號(hào)與噪聲疊加。對(duì)加噪后的信號(hào)做出時(shí)域波形和頻譜圖，觀察波形的變化。這里要求設(shè)計(jì)數(shù)字濾波器，根據(jù)對(duì)加噪前的頻譜以及加噪后的頻譜的觀察，可以采用低通考慮設(shè)計(jì)相對(duì)較為簡(jiǎn)單的巴特沃斯低通濾波器進(jìn)行濾波。濾波前首先要確定設(shè)計(jì)指標(biāo)，觀察頻譜這里暫取0然后開始設(shè)計(jì)巴特沃斯低通濾波器,這里我把設(shè)計(jì)的程序打

11、包成一個(gè)函數(shù)方便調(diào)用，函數(shù)的框圖如圖三（巴特沃斯低通濾波器開始k= 1 , n參量輸入qc = ftp(10°lap-1)nc = np(io°lap 由rz,pz,kz合幷部分分式得hr(z)由rz,pz,kz合幷部分分式得hr(z)2avs(k) = nrp(fc)jf«(s項(xiàng)式合并確定九分子多項(xiàng)式民(s)部分分式展開得rs.psksz(k)=訂只如rss=rs/fsx;>pz(k.l)=iyz(k)nlength b < n+1| b(k) = 0(*!=；bz =bz= mzla| az=flip(a)n輸出占乙/結(jié)束圖三 巴特沃斯低通濾波器函

12、數(shù)運(yùn)算流圖)可以發(fā)現(xiàn)fft的運(yùn)算過程很有規(guī)律性。h二2"點(diǎn)的m共進(jìn)行m級(jí)運(yùn)算，每級(jí)由n/2個(gè)a(0)a(0)a(0)a(0)x(0)x(1)xxxxxx8點(diǎn)dit-fft運(yùn)算流圖圖四函數(shù))所示。該函數(shù)調(diào)用需要5個(gè)輸入?yún)?shù)分別為：n.：模擬通帶邊界角頻率，qs：模擬 阻帶邊界角頻率，耳：通帶最大衰減，l：阻帶最小衰減，fs：信號(hào)釆樣頻率。參數(shù)輸入后，根據(jù)數(shù)學(xué)上的推導(dǎo)開始計(jì)算。首先根據(jù)書本公式6.2.18確定濾波器的階數(shù) n,然后根據(jù)書本公式6213求出歸一化極點(diǎn)卩卜 接著求出3db截止頻率qy然后對(duì)卩£進(jìn)行去歸一化 得到» 然后得到h/s)的分子分母系數(shù)，接著對(duì)進(jìn)行

13、部分分式展開，由書本公式6.3.6, 利用脈沖響應(yīng)不變法將f4(s)變換到h(z)，接著合并部分分式，求得 h(z)的分子分母多項(xiàng)式的系數(shù)，對(duì)多項(xiàng)式的系數(shù)進(jìn)行整理排列得到巴特沃斯低通濾波器系統(tǒng)函數(shù)的分子分母多項(xiàng)式的系數(shù)向量b乙az,作為函數(shù)輸出。接著調(diào)用freqz函數(shù)求出設(shè)計(jì)的低通濾波器的頻率 響應(yīng)特性函數(shù)，然后繪制出去頻響特性曲線并觀察，看是否符合設(shè)計(jì)指標(biāo)。若符合則開始 濾波，不符合就要重新設(shè)計(jì)。濾波時(shí)直接調(diào)用了 filter函數(shù)，得到的信號(hào)做出時(shí)域波形和頻 譜圖并與濾波前的信號(hào)對(duì)比，觀察濾波效果，是否濾除了100khz的干擾信號(hào)，并保留了原有的音頻信號(hào)。實(shí)驗(yàn)內(nèi)容三：編寫 fft算法程序此部

14、分要求自己編寫fft算法程序，首先要充分了解fft的運(yùn)算過程。觀察圖四(8點(diǎn)dit-fft蝶形運(yùn)算組成。同一級(jí)中每個(gè)蝶形的兩個(gè)輸入數(shù)據(jù)只對(duì)計(jì)算本蝶形有用，而且每個(gè)蝶形的輸入、輸出數(shù)據(jù)結(jié)點(diǎn)又同在一條水平線上，這意味著計(jì)算完一個(gè)蝶形后，所得輸出數(shù)據(jù)可立即存入原輸入數(shù)據(jù)所占用的存儲(chǔ)單元。這樣，經(jīng)過 m級(jí)運(yùn)算后，原來存放輸入序列數(shù)據(jù)的n個(gè)存儲(chǔ)單元(數(shù)組 a)中便一次存放 x(k)的n個(gè)值。這樣可以節(jié)省大量?jī)?nèi)存，從而降低設(shè)備成本。轉(zhuǎn)因子的指數(shù)。但各級(jí)的旋轉(zhuǎn)因子和循環(huán)方式都有所不同。為了編寫計(jì)算程序，應(yīng)先找出旋轉(zhuǎn)因子與運(yùn)算級(jí)數(shù)的關(guān)系。用表示從左到右的運(yùn)算級(jí)數(shù)。 不難發(fā)現(xiàn)，第l級(jí)共有p :個(gè)不同的旋轉(zhuǎn)因子。

15、對(duì)一般情讀入x(n),m況，第 l級(jí)的旋轉(zhuǎn)因子為昭遡斛冃仏嚴(yán)因?yàn)?、山尸二憶曲所這樣就確定了第 l級(jí)運(yùn)算的旋轉(zhuǎn)因子。接著研究蝶形運(yùn)算規(guī)律，設(shè)序列 x(n)經(jīng)吋域抽選(倒序)后，按圖中次序存入數(shù)組a中，如果蝶形運(yùn)算的兩個(gè)輸入數(shù)據(jù)相距b個(gè)點(diǎn)，應(yīng)用原位計(jì)算，則蝶形運(yùn)算相距b個(gè)點(diǎn)，應(yīng)用原位計(jì)算，則蝶形運(yùn)算開始倒序j= o.b-1表示成如下形式心(/)"“(/)+也/+3)昭 al(j+b)t = a(k)+a(lc + b)昭/ = 0x2.m,2w-1；下標(biāo)l表示第l級(jí)運(yùn)算，4(/)則表示第l級(jí)運(yùn)算后的數(shù)組元素 a(j)的值。ffu-10)表示第l輸出結(jié)束級(jí)運(yùn)算前a(j)的值。圖五 fft

16、運(yùn)算程序框圖另外由圖四還有一些運(yùn)算規(guī)律：第l級(jí)中，每個(gè)蝶形的兩個(gè)輸入數(shù)據(jù)相距3二?個(gè)點(diǎn)；每級(jí)有b個(gè)不同的旋轉(zhuǎn)因子；同一旋轉(zhuǎn)因子對(duì)應(yīng)著間隔為卩點(diǎn)的個(gè)蝶形?？偨Y(jié)上述規(guī)律，采用以下算法。先從輸入端開始，逐級(jí)進(jìn)行，共進(jìn)行m級(jí)運(yùn)算。在進(jìn)行第 l級(jí)運(yùn)算時(shí)，依次求岀 b個(gè)不個(gè)蝶形。這樣，利用三重循環(huán)實(shí)同的旋轉(zhuǎn)因子，每求出一個(gè)旋轉(zhuǎn)因子，就計(jì)算完它對(duì)應(yīng)的所有 現(xiàn)fft算法，程序框圖如圖五(fft運(yùn)算程序框圖)所示。忙外，ditffti云算流圖的輸出 x(k)為白然順丿命，但一為了i舌應(yīng)原付計(jì)算，其輸入丿命列不是按的自然順序排列，憩建m次偶奇抽選后的排序稱為序列x(n)的倒序。因此，在運(yùn)算 m級(jí)蝶形之前應(yīng)對(duì)序列x

17、(n)進(jìn)行倒序。由于本程序使用的是 matlab高級(jí)語言無法直接倒置二進(jìn)制數(shù)位， 因 此必須找到產(chǎn)生倒序的十進(jìn)制運(yùn)算規(guī)律。觀察倒序的二進(jìn)制數(shù)可以發(fā)現(xiàn)，倒序數(shù)是在m位二進(jìn)制數(shù)最高位加1,逢2向低位進(jìn)位。為了敘述方便，用j表示當(dāng)前倒序數(shù)的十進(jìn)制數(shù)值。對(duì)于h二2" m位二進(jìn)制數(shù)最高位的十進(jìn)制權(quán)值為n/2,lh=n/2j =lh2 =n-21且從左向右二進(jìn)制位的權(quán)值依次為n/4, n/8,，2, 1o因此，最高位加1相當(dāng)于十進(jìn)制運(yùn)算 j+n/2。如果最高位是0 (jvn/2),則直接由j+n/2得下一個(gè)倒序值；如果最高位是1 ( j> n/2),則先將最高位變成 。問7/2),i &g

18、t; j ?ni = 1 , n1然后次高位加 1(j+n/4) o但次高位加1時(shí)，同樣要斷0、1值；依次類推直到完成最高位加逢2向右進(jìn)位的運(yùn)算。倒序的程序框圖如圖六 (倒序程序框圖)所示。形成倒序j后，將元數(shù)組中存放的輸入序列重新按倒序排序。 設(shè)原輸入序列x(l)先按自然順序存放入數(shù)組n?j = j-kk = k/2t = a(l) a(l) = a(j)進(jìn)行倒序，發(fā)現(xiàn)第一個(gè)勺最后一個(gè)虛列值不需要重排， 所以順序數(shù)圖六倒序程序框圖中。對(duì)x(n)i的初值為1,盤為n-2,倒序數(shù)j的初值為n/2o每計(jì)算出一個(gè)倒序值j便與循環(huán)語句自動(dòng)生成的順序i比較,當(dāng)i二j吋，不需要交換，當(dāng) 岡吋，a和a(j)

19、交換數(shù)據(jù)。另外，為了避免重復(fù)調(diào)換，對(duì)ivj的情況進(jìn)行調(diào)換。以上是fft函數(shù)的編程思路，具體程序見附錄，已純yj_fft(函數(shù)，方便調(diào)用，附錄中遞對(duì)同一個(gè)音頻序列用 yj_fft(和fft()進(jìn)行傅里葉變換的比較測(cè)試程序，用于觀察設(shè)計(jì)的函數(shù)的數(shù) 性。八、實(shí)驗(yàn)任務(wù)windows系統(tǒng)中的“ding喑頻信號(hào)的采集、分析、合成(程序框圖見圖一）contents.音頻信號(hào)的采集12音頻信號(hào)的頻譜分析自動(dòng)搜尋頻譜的峰值點(diǎn)并在圖中表示出來1.3.音頻信號(hào)的分解和合成1.4音頻信號(hào)的回放1.5音頻信號(hào)分段傅里葉分析% % title:數(shù)字信號(hào)處理實(shí)驗(yàn)一一音頻信號(hào)采集、分析及處理% author : 120421

20、101 楊杰% % % 實(shí)驗(yàn)任務(wù)一：windows系統(tǒng)中的“ding”音頻信號(hào)的采集、分析、合成% %clear all ;close all ;1.1 音頻信號(hào)的采集sign al,fs = audioread('ding.wav* );淤取音頻文件，得到音頻序列signal及采樣頻率fssignal = sign al(:,1);%只取第一個(gè)聲道的序列l(wèi)en_signal = length signal);% 言號(hào)的長(zhǎng)度t= 1/fs;嫁樣周期t = t:t:le n_ signal/fs;%!寸間t序列sound( signal,fs);%原音頻信號(hào)播放pause(1);痼停1s

21、1.2.音頻信號(hào)的頻譜分析figure 'name*音頻信號(hào)的頻譜分析'position ,50,350,400,300）;涯成一個(gè)窗口subplot( 2,1,1 );劉各窗口分為兩行，第一部分用于繪制時(shí)域波形plot( t,signal);滋制時(shí)域波形axis( t len_signal/fs -0.08 0.08);%艮定顯示范圍title（|時(shí)域波形圖'）;淞制該波形的標(biāo)題xlabel( 1 時(shí)間(s); ylabel(1 幅值')；滋制橫縱坐標(biāo)的標(biāo)簽fft_n = 2anextpow2( len_signal);加十算fft變換的點(diǎn)數(shù)signal =

22、fft( signal,ffi_n );%快速傅里葉變換signal.amp = abs( signal );%只做幅頻分析，對(duì)fft變換后序列取模f = fs*(0:fft_n/2-1)/fft_n;瀰率橫軸f,由于對(duì)稱性，只取前一半signal amp h/lfv4gn/l amp (斑%subplot( 2,1,2 );%窗口第二部分，繪制頻譜圖08唸會(huì)制頻譜圖唏示題囁 會(huì)制橫縱坐標(biāo)的標(biāo)簽?zāi)聂薅@示范圍plot( f,signal_amp_half );title('頻譜圖)；xlabel( * 頻率(hz );ylabel( * 幅度); axis( 0 4000 0 250)

23、;時(shí)域波形圖0 050-0.050406時(shí)間 頻譜圖b頻率(hz)下面程序?qū)π盘?hào)的分解分析：觀察音頻信號(hào)的頻譜圖可知該音頻信號(hào)能量主要集中在兩個(gè)頻率上, 主要是找到這兩個(gè)峰值點(diǎn)對(duì)應(yīng)的頻率。自動(dòng)搜尋頻譜的峰值點(diǎn)并在圖中表示出來signal.pksjocs = findpeaks( signal_amp_half ); sort_pks = sort( signal_pks, 'descend'); sortjoc = find( signal_pks>=sort_pks(2);rd 二 locs( sort_loc(1);n2 = locs( sort_loc(2);y1

24、= signal_amp_half(n1);y2 = signal_amp_half(n2);f1 = (n1-1)*fs/fft_n;f2 = (n2-1)*fs/fft_n;s1 = sprintf(* 峰值一 (%fhz)' ,f1 );s2 = sprintf(* 峰值二(%fhz)' ,f2);text( f1,y 1,s1 );text( f2,y2,s2 );%尋找頻譜峰值點(diǎn)加寸頻譜峰值點(diǎn)進(jìn)行排序呀戈到峰值最大的2個(gè)點(diǎn)排序后在pks中的位置艱高峰值點(diǎn)在序列中的位置癖二高峰值點(diǎn)在序列中的位置艱高峰值點(diǎn)在頻譜中的幅度瞬二高峰值點(diǎn)在頻譜中的幅度鴻高峰值點(diǎn)頻率獅二高峰值點(diǎn)

25、頻率漩義字符串1漩義字符串2陸圖上最高峰值點(diǎn)處顯示“峰值一”拓圖上第二高峰值點(diǎn)處顯示“峰值二”時(shí)域波形圖0.05 0-0.05020406時(shí)間0.8頻率(hz)分析：找到了兩個(gè)峰值點(diǎn)對(duì)應(yīng)的頻率并在圖上顯示出來。1.3.音頻信號(hào)的分解和合成y1 =(y1/(y1+ y2)*cos(2*pi*f1*t); y2 = (y2/(y1+ y2)*cos(2*pi*f2*t);瞬一個(gè)頻率信號(hào)瞬二個(gè)頻率信號(hào)a = max( signal);hecheng = (y1 +y2)/max( y1+y2 );hecheng = hecheng*a;%原信號(hào)的最大幅度殆成信號(hào)歸一化酬吏合成信號(hào)與原信號(hào)的幅度相等f

26、igure( name*音頻信號(hào)的合成，position1 subplot( 2,1,2 );plot( t,hecheng );title( *合成音頻時(shí)域波形圖)；xlabel( * 時(shí)間(s) ); ylabel，幅值)；,450,350,400,300);儺成一個(gè)窗口加各窗口分為兩行，第二部分用于繪制時(shí)域波形呀會(huì)制合成信號(hào)時(shí)域波形唏會(huì)制該波形的標(biāo)題淤制橫縱坐標(biāo)的標(biāo)簽subplot( 2,1,1 );h1=plot( t,signal);axis( t len_signal*t -0.08 0.08); title('原信號(hào)時(shí)域波形圖1 );xlabel( * 時(shí)間(s) );

27、ylabel( 1 幅值);煖口第一部分，繪制原信號(hào)時(shí)域波形喏會(huì)制原信號(hào)時(shí)域波形哪艮定顯示范圍滋制標(biāo)題淞制橫縱坐標(biāo)的標(biāo)簽figure( 'name* ,'合成音頻時(shí)域波形圖(時(shí)域放大) plot( t,hecheng );title('合成音頻時(shí)域波形圖(時(shí)域放大)；xlabel('時(shí)間(s)'); ylabel( '幅值)；axis 0.2 0.21 -0.08 0.08);1 , position* ,850,350,400,300);成一個(gè)窗口淞制時(shí)域波形淞制標(biāo)題淞制橫縱坐標(biāo)的標(biāo)簽滋艮定顯示范圍(局部放大)%生成一個(gè)窗口figure( &#

28、39;name'， 音頻信號(hào)的分解,'position1 ,50,10,400,300);subplot( 2,1,1 );fft_y1 = abs( fft( y1*a,fft_n );plot( f,fft_y11:(fft_n/2);title('第一個(gè)頻譜峰值點(diǎn)信號(hào)頻率譜)；xlabel( * 頻率(hz);text( f1,fft_y1(n1),s1 );axis 0 4000 0 1000);subplot( 2,1,2 );fft_y2 = abs( fft( y2*a,fft_n );plot( f,fft_y2(1:(fft_n/2);title(&#

29、39;第二個(gè)頻譜峰值點(diǎn)信號(hào)頻率譜)；xlabel( * 頻率(hz );ylabel( * 幅度);text( f2,fft_y2(n2),s2);axis( 0 4000 0 1000);figure( name*合成前后頻率譜比較，position* subplot( 2,1,2 );hecheng = abs( fft(hecheng,fft_n);plot( f,hecheng(1:(fft_n/2);title('合成后信號(hào)頻率譜-)；xlabel( 1 頻率(hz );ylabel'幅度)；text( f1,fft_y1(n1),s1 );text( f2,fft_

30、y2(n2),s2);axis 0 4000 0 1000);subplot( 2,1,1 );plot( f,signal_amp(1:(fft_n/2);title('原信號(hào)頻譜1 );xlabel('頻率(hz );ylabel'幅度)；axis 0 4000 0 250);text( f1,y1,s1 );text( f2,y2,s2 );剌各窗口分為兩行，第一部分用于繪制第一個(gè)頻率信號(hào)的頻譜%對(duì)第一個(gè)頻率信號(hào)進(jìn)行 fft ,并取模 您制頻譜圖%繪制標(biāo)題淹制橫縱坐標(biāo)的標(biāo)簽%在峰值點(diǎn)處顯示“峰值一” 卅艮定顯示范圍爛口第二部分，第二個(gè)頻率信號(hào)的頻譜 慟第二個(gè)頻率信

31、號(hào)進(jìn)行fft ,并取模 淞制頻譜圖%會(huì)制標(biāo)題滋制橫縱坐標(biāo)的標(biāo)簽獨(dú)峰值點(diǎn)處顯示“峰值二” 哪艮定顯示范圍,450,10,400,300);淞成一個(gè)窗口爛口第二部分，繪制合成信號(hào)頻譜圖 腋寸合成信號(hào)進(jìn)行fft變換并取模 燈制頻譜圖淞制標(biāo)題濰制橫縱坐標(biāo)的標(biāo)簽峰值點(diǎn)處顯示“峰值一” 憾峰值點(diǎn)處顯示“峰值二” 滋艮定顯示范圍嘀口第一部分，繪制原信號(hào)頻譜圖 淞制原信號(hào)頻譜圖滋制標(biāo)題滋制橫縱坐標(biāo)的標(biāo)簽%艮定顯示范圍濰圖上最高峰值點(diǎn)處顯示“峰值一” 注圖上第二高峰值點(diǎn)處顯示“峰值二”原信號(hào)時(shí)域波形圖時(shí)間合成音頻時(shí)域波形圖010.2040.60.8時(shí)間.20 2020 2040 2060 208時(shí)間合成音頻時(shí)域

32、波形圖（時(shí)域放大）0.080.060.040.02凹 00.02-0.04-0.060峠分析：觀察音頻信號(hào)合成后的時(shí)域波形可以發(fā)現(xiàn)合成后的音頻信號(hào)與原始音頻信號(hào)差距較大，這是因?yàn)樵家纛l信號(hào)從頻譜圖上看，雖然看似能量都集中在兩個(gè)頻率上，但實(shí)際上，在那兩個(gè)頻率的 兩旁還有許多頻率分量，用這種音頻信號(hào)合成方法合成音頻信號(hào)時(shí)，忽略掉了其他的所有頻率分量， 合成出來的音頻信號(hào)自然與原信號(hào)有很大差異。第一個(gè)頻謚峰值點(diǎn)信號(hào)頻率譜顎率(hz)第二個(gè)頻語峰值點(diǎn)信號(hào)頻率譜1000 黒 500梅值二(1242 869568hz)-90 1 01000200030004000頻率(hz)原信號(hào)頻譜200磬100峰値

33、一(93131035hz) 殍值二(1242 869568hz)°01000100020003000頻率(hz)合成后信號(hào)頻率譜予伯35甩)_l4000值二(1242 869568hz)°0100020003000頻率(hz)4000%原音頻信號(hào)播放瞬停1s殆成音頻信號(hào)播放瞬停1s可以發(fā)現(xiàn)分析這是分解的音頻信號(hào)的頻譜圖以及合成后的頻譜圖和原信號(hào)的頻譜圖的比較。計(jì)算得到的兩個(gè)頻譜峰值對(duì)應(yīng)的頻率是正確的。1.4音頻信號(hào)的回放sound( signal,fs); pause(1);sound( hecheng,fs); pause(1);1.5音頻信號(hào)分段傅里葉分析n = 102

34、4;m = ceil( len_signal/n );fftfenduan 二 zeros( m,n );for k =fftfenduanf k,:) = fft( siqnal(n*(k-1)+1:n*k),n );滋定每段的點(diǎn)數(shù)2切，此處為1024漏定分的段數(shù)炬成n行m列數(shù)組用于存儲(chǔ)分段fft后的序列%分段傅里葉變換endfftfenduan( m,:) = fft( signal(n*k+1:len_signal),n );fenduan = zeros( 1,m*n );for k = 1:m;fenduan(1,n*(k-1)+1:n*k) = ifft( fftfenduan(k

35、,:),n );end演后一行音頻信號(hào)進(jìn)行 m%生成m5l長(zhǎng)序列，用于存儲(chǔ)ifft 后的序列%分段傅里葉逆變換并合成%分段傅里葉分析后的音頻信號(hào)的播放,850,10,400,300);炬成一個(gè)窗口頰口第一部分，繪制分段傅里葉分析后信號(hào)時(shí)域波形 %會(huì)制分段傅里葉分析后信號(hào)時(shí)域波形哪艮定顯示范圍強(qiáng)會(huì)制標(biāo)題唸會(huì)制橫縱坐標(biāo)的標(biāo)簽滋寸分段傅里葉分析后信號(hào)做fft ,并取模煖口第二部分，繪制分段傅里葉分析后信號(hào)頻譜圖強(qiáng)會(huì)制分段傅里葉分析后信號(hào)頻譜圖唸會(huì)制標(biāo)題您制橫縱坐標(biāo)的標(biāo)簽?zāi)聂薅@示范圍sound( fenduan,fs );figure( name，,'分段傅里葉分析后時(shí)域頻域圖，positi

36、on*subplot(2,1,1);plot( t,fenduan(1:len_signal);axis t len_signal/fs -0.08 0.08);title(*分段傅里葉分析后音頻時(shí)域波形圖)；xlabel( * 時(shí)間(s); ylabel( * 幅值)；%分段分析合成后的傅里葉變換頻譜fenduan = abs( fft( fenduan,fft_n );subplot(2,1,2);plot( f,fenduan(1:floor(fft_n/2);title('分段分析后音頻頻率譜* );xlabel('頻率(hz );ylabel'幅度)；axis

37、( 0 4000 0 250);分段傅里葉分析后音頻時(shí)域波形圖0 050-0.050.40.6時(shí)間08頻率(hz)分析：這是分段傅里葉分析后合成的音頻信號(hào)的時(shí)域波形及頻譜圖，發(fā)現(xiàn)和原始信號(hào)一樣，相對(duì)于 一次性傅里葉分析，很好地實(shí)現(xiàn)了信號(hào)的合成。因?yàn)橐淮涡愿道锶~分析過程中只考慮了單一的兩個(gè) 頻率分量，忽略了其他所有的頻率分量，而分段傅里葉分析考慮了所有的頻率?？偨Y(jié)，對(duì)于音樂或 說話的聲音，一般不可能只有單一的幾個(gè)頻率合成，它們往往是由許許多多的頻率分量組成的，考 慮的頻率越多，所合成的信號(hào)越接近原始的音頻信號(hào)。published with matlab? r2014a實(shí)驗(yàn)內(nèi)容二：任意音頻信號(hào)的時(shí)

38、域和頻域分析及數(shù)字濾波器設(shè)計(jì)（程序框 圖見圖二）contents2.1 音頻信號(hào)的采集22音頻信號(hào)的頻譜分析23引入干擾噪聲100k2.4.數(shù)字濾波器設(shè)計(jì)25用數(shù)字濾波器進(jìn)行濾波調(diào)用flite函數(shù)% title :數(shù)字信號(hào)處理實(shí)驗(yàn)一一音頻信號(hào)采集、分析及處理% author : 120421101 楊杰 % % % 實(shí)驗(yàn)任務(wù)二：任意音頻信號(hào)的時(shí)域和頻域分析及數(shù)字濾波器設(shè)計(jì)% %clear all ;close all ;2.1 音頻信號(hào)的采集sign al,fs = audioread('reny i.wav1);%讀取音頻文件，得到音頻序列signal及采樣頻率fslen_sig n

39、al = fs;signal = signal(1:le resign al,1);%截取信號(hào)的長(zhǎng)度為1s%只取第一個(gè)聲道的序列t= 1/fs;殊樣周期t = t:t:len_s ign al*t;%時(shí)間序列t,50,350,400,300);2.2.音頻信號(hào)的頻譜分析figure( 'name'， 音頻信號(hào)分析','position* subplot( 2,1,1 );plot( t,signal);title('時(shí)域波形圖* );xlabel( 1 時(shí)間,);ylabel( 1 幅值')；fft_n = 2anextpow2( len_sig

40、nal);signal = fft( signal,fft_n );signal.amp = abs( signal );f = fs*(0:fft_n/2-1)/fft_n;signal_half 二 signal_amp(1:(fft_n/2);%生成一個(gè)窗口%各窗口分為兩行，第一部分用于繪制時(shí)域波形%繪制時(shí)域波形%繪制該波形的標(biāo)題%繪制橫縱坐標(biāo)的標(biāo)簽硏算fft變換的點(diǎn)數(shù)漱速傅里葉變換頫做幅頻分析，對(duì)ffi變換后序列取模獅率橫軸f,由于對(duì)稱性，只取前一半%由于對(duì)稱性丿貝取前一半subplot( 2,1,2 );滴口第二部分，繪制頻譜圖嗨會(huì)制頻譜圖術(shù)題滋制橫縱坐標(biāo)的標(biāo)簽%艮定顯示范圍%提高信

41、號(hào)的采樣頻率10倍，以滿足可以引入fs = fs*10;tt= 1/fs;tt = t:tt:len_signal*t;signalfs = interpl (t,signal,tt); len_signalfs = length( signalfs );figure( name，,'提高采樣頻率后的音頻信號(hào) subplot( 2,1,1 );plot( tt,signalfs);title(1提高采樣頻率后音頻時(shí)域波形圖xlabel( * 時(shí)間(s); ylabel( * 幅值)；100khz的噪聲%各采樣頻率提高10倍 淞高采樣頻率后的采樣周期 加寸間序列tt瀏原信號(hào)一階線型插值

42、漩值后的音頻信號(hào)的長(zhǎng)度 , position* ,450,350,400,300);池成一個(gè)窗口剌各窗口分為兩行，第一部分用于繪制采樣頻率提高后的時(shí)域波形滋制釆樣頻率提高后的時(shí)域波形嚇題強(qiáng)會(huì)制橫縱坐標(biāo)的標(biāo)簽fft_nn = 2anextpow2( len_signalfs );signalfs = fft( signalfs,fft_nn );signalfs.amp = abs( signalfs );ff = fs*(0:fft_nn/2-1)/fft_nn;signalfs_amp_half = signalfs_amp( 1:(fft_nn/2);subplot( 2,1,2 );pl

43、ot( ff,signalfs_amp_half );title('提高采樣頻率后音頻頻譜)；xlabel( * 頻率(hz );ylabel( * 幅度);axis( 0 5000 0 4000);力十算采樣頻率提高后fft變換的點(diǎn)數(shù) 於夬速傅里葉變換頫做幅頻分析，對(duì)fft變換后序列取模 獅率橫軸ff ,由于對(duì)稱性，只取前一半 %由于對(duì)稱性，只取前一半癖二部分用于繪制采樣頻率提高后的頻譜圖滋制采樣頻率提高后的頻譜圖嚇題淞制橫縱坐標(biāo)的標(biāo)簽加艮定顯示范圍0020.40.6081時(shí)間頻譜plot( f,signal_half );title('頻譜)；xlabel( * 頻率(hz

44、 );ylabel( * 幅度); axis( 0 5000 0 400);頻率(hz)0.2提高采樣頻率后音頻時(shí)域汲形團(tuán)-0.20020.406081時(shí)間提高采祥頻率后音頻頻諳劭率(hz)分析：首先看原信號(hào)的時(shí)域波形及頻譜圖，可知該信號(hào)截取的是 1s,頻率分量主要集中在低頻段, 提高采樣頻率后比較時(shí)域波形及頻譜圖還是和原始信號(hào)一樣，說明提高采樣頻率過程中并沒有改變 信號(hào)，只是提高了采樣頻率，符合要求。2.3.引入干擾噪聲100knoise = 0.5*cos(2*pi*100*10a3*tt);signal_noise = sign alfs + no ise1:le ngth(signal

45、fs); figure( 'name* , 1引入噪聲后的音頻信號(hào),'position* subplot( 2,1,1 );plot( tt,signal_noise);title('引入噪聲后音頻時(shí)域波形圖)；xlabel( * 時(shí)間(s) ); ylabel，幅值)；%100khz 噪聲啕入噪聲,850,350,400,300);淞成一個(gè)窗口劉各窗口分為兩行，第一部分用于繪制引入噪聲后的時(shí)域波形 啕【入噪聲后的時(shí)域波形 滋示題淤制橫縱坐標(biāo)的標(biāo)簽加寸加噪后信號(hào)進(jìn)行fft頫做幅頻分析，對(duì)fft 變換后序列取模%由于對(duì)稱性，只取前一半瞬二部分用于繪制加噪聲后音頻頻率譜 淞

46、制加噪聲后音頻頻率譜 術(shù)題滋制橫縱坐標(biāo)的標(biāo)簽signal.noise = fft( signal_noise,fft_nn );signal_noise_amp = abs( signal.noise );signalfs_noise_amp_half = signal_noise_amp( 1:(fft_nn/2); subplot( 2,1,2 );plot( ff,signalfs_noise_amp_half );title('加噪聲后音頻頻譜1 );xlabel( * 頻率(hz)，);ylabel( * 幅度);引入哩聲后音頻時(shí)域波形圉liiii"1o021020

47、406時(shí)間加噪聲后音頻傾譜0.0°0468頻率(hz)10 12x1owp = 2*pi*4000;ws = 2*pi*8000;ap = 1；as = 60;b乙az = yj_butter_low(wp,ws,ap,as,fs);h,w = freqz( b乙az);hf = abs( h );figure( name，,'巴特沃斯低通濾波器系統(tǒng)函數(shù)圖 subplot( 2,1,1 );plot( w,hf);title('濾波系統(tǒng)幅頻特性曲線* );xlabel('數(shù)字頻率(0兀);ylabel( * 幅度,); subplot( 2,1,2);plot

48、( w,20*log10(hf);title('濾波系統(tǒng)幅頻特性曲線(db) * ); xlabel(數(shù)字頻率(0龍)*); ylabel( 'db');滅離散系統(tǒng)頻率響應(yīng)特性函數(shù) 頒模分析：這是加噪后的信號(hào)時(shí)域波形及頻譜圖，可以知道該信號(hào)時(shí)域波形發(fā)生了較大的改變，從頻譜 圖中也可以發(fā)現(xiàn)在10人5處即1ookhz處出現(xiàn)了一個(gè)較高的沖激。說明成功引入了 1ookhz的正弦干 擾信號(hào)。2.4.數(shù)字濾波器設(shè)計(jì)%巴特沃斯低通濾波器帶邊界角頻率滋且?guī)н吔缃穷l率帶最大衰減% 且?guī)ё畲笏p瀾自己設(shè)計(jì)的巴特沃斯低通濾波器函數(shù)求系統(tǒng)函數(shù)分子分母 , position* ,50,10,40

49、0,300）;%成一個(gè)窗口滋各窗口分為兩行，第一部分用于繪制濾波系統(tǒng)幅頻特性曲線滋會(huì)制濾波系統(tǒng)幅頻特性曲線 悄示題滋制橫縱坐標(biāo)的標(biāo)簽癖二部分用于繪制濾波系統(tǒng)幅頻特性曲線（db）燈制濾波系統(tǒng)幅頻特性曲線（db）漪題淤制橫縱坐標(biāo)的標(biāo)簽遞波系統(tǒng)幅頻特性曲線000.511.522533.5數(shù)字頻率(hrr數(shù)字頻率orr分析：這是設(shè)計(jì)的巴特沃斯低通濾波器的幅頻特性曲線，觀察該曲線可知該系統(tǒng)屬于低通濾波器, 且符合設(shè)計(jì)指標(biāo)的要求。2.5用數(shù)字濾波器進(jìn)行濾波 調(diào)用filter函數(shù)溺用函數(shù)濾波sign aljvbo = filter(b z,az,sig nal_noise);figure( name*，濾波前后時(shí)域信號(hào)比較，卩osition* subplot(2,1,2);plot(tt,signal_lvbo);title('濾波后音頻波形1 );xlabel( * 時(shí)間(s)' );ylabel(* 幅值')；subplot(2,1,1);plot( tt,signal_noise);title('濾