現(xiàn)代仿真作業(yè)

1、MATLAB在音頻信號處理中的應(yīng)用1.引言 MATLAB是美國Math Works公司推出的一種面向工程和科學(xué)計算的交互式計算軟件、它以矩陣運算為基礎(chǔ)，把計算、可視化、程序設(shè)計融合在一個簡單易用的交互式工作環(huán)境中，是一款數(shù)據(jù)分析和處理功能都非常強大的工程實用軟件。MATLAB軟件除可以對文本數(shù)據(jù)進行操作外，還可以對圖像文件、WAV類型的音頻文件和AVI類型的視頻文件進行讀寫.借助windows系統(tǒng)自帶的音頻錄制播放設(shè)備和MATLAB軟件，可以完成音頻信號的綜合分析處理。音頻是信號的一種，處理數(shù)字音頻信號也是一種數(shù)字信號分析與處理。人能聽到的聲音頻率范圍為20Hz到20kHz，語音信號頻率為30

2、0Hz到3.4kH。木文首先對音頻信號的采集與頻譜分析作介紹，然后介紹了用MATLAB處理音頻信號的整體流程，并且對我們平常所聽的音樂文件進行一系列的處理，然后用我們的耳朵親自去聽，去感受不同處理方法的效果。2.音頻信號的采集與頻譜分析進行音頻的頻譜分析時，首先要對聲音信號進行采。MATLAB 的數(shù)據(jù)采集工具箱提供了一整套命令和函數(shù),通過調(diào)用這些函數(shù)和命令,可直接控制聲卡進行數(shù)據(jù)采集。Windows 自帶的錄音機程序也可驅(qū)動聲卡來采集語音信號，并能保存為WAV 格式文件供MATLAB 相關(guān)函數(shù)直接讀取、寫入或播放。函數(shù)y, fs, bits =wavre

3、ad(Blip, N)。用于讀取音頻，對音頻信號進行采集，采樣值放在參數(shù)y中，fs表示每秒采樣點數(shù)，即采樣頻率，bits表示每個采樣點在編碼時所占位數(shù)。N表示采樣點總數(shù)。參數(shù)Blip為音頻所在地址，如:'C:yinpinl'。調(diào)用函數(shù)fft可對己采集音頻信號進行時頻轉(zhuǎn)換，通過函數(shù)abs()和angle()可分別得到信號頻譜的幅頻圖和相頻圖。對放在C盤目錄下手機短信鈴聲”test.wav”進行采樣與頻譜分析。其體代碼如下:clear all; y0,Fs0,nbits0=wavread('test.wav'); a=round(length(y0); y01=ff

4、t(y0);y011=abs(y01); t0=linspace(0,a/Fs0,a); figure(1);subplot(2,1,1);plot(t0,y0); title('(a)','fonts',10.5 ,'position' ,25,-1.5); xlabel('時間/s','fonts',10.5 ,'position',45,-1.4); ylabel('幅值n','fonts', 10.5,'position',-3,0.7); f

5、0=linspace(0,Fs0,a); figure(2);subplot(2,1,1);plot(f0,y011); title('(a)','fonts',10.5 ,'position' ,8000,-780); xlabel('頻率f/hz','fonts',10.5 ,'position',14000,-600); ylabel('幅值n','fonts', 10.5,'position',-1500,2000); figure(3);su

6、bplot(2,1,1); plot(f0(1:round(a/16),y011(1:round(a/16); title('(a)','fonts',10.5 ,'position' ,500,-780); xlabel('頻率f/hz','fonts',10.5 ,'position',850,-600); ylabel('幅值n','fonts', 10.5,'position',-90,2000); 運行程序，得到音頻信y0，采樣率Fs0 =4

7、4100/s;采樣總點數(shù)a= 441000;音頻總時長t0=10s。通過圖1音頻信號時域波形圖，可以看出音頻頻率的大小，但不夠定量，頻率分布定量計算需要通過頻譜分析。通過MATLAB對音頻進行處理，得到圖1的音頻信號時域信號波形圖，和圖2、圖3的幅頻圖。圖1音頻信號時域信號波形圖圖2音頻信號幅頻圖圖3音頻信號頻譜密集部分幅頻圖由圖2的音頻信號幅頻圖可以看圖，此音頻信號的頻率主要集中在低頻區(qū)和高頻區(qū)，此外我們還可以用MATLAB計算頻率在1000hz以下的頻譜部分占0-45000hz頻譜部分的比例。具體實現(xiàn)語句為：w0=(sum(y011(1:10000)/sum(y011(1:45000);運

8、行結(jié)果w0=0.5324。3.MATLAB 處理音頻信號的流程我們聽到的聲音一般都屬于復(fù)音,其聲音信號由不同的振幅與頻率的波合成而得。WAV 格式音頻信號作為分析處理的輸入數(shù)據(jù)。用MATLAB 處理音頻信號的基本流程是先將WAV 格式音頻信號經(jīng)wavread 函數(shù)轉(zhuǎn)換成MATLAB 列數(shù)組變量，再用 MATLAB 對音頻信號進行音量標準化，聲道分離合并與組合，數(shù)字濾波，數(shù)據(jù)轉(zhuǎn)換 。處理后的數(shù)據(jù)如是音頻數(shù)據(jù)則可用wavwrite 轉(zhuǎn)換成 WAV 格式文件或用sound wavplay&#

9、160;等函數(shù)直接回放。3.1音量標準化錄制聲音過程中需對聲音電平進行量化處理，最理想的量化是最大電平對應(yīng)最高量化比特。但實際卻很難做到，常有音輕問題。利用MATLAB很容易實現(xiàn)音量標準化。即最大電平對應(yīng)最高量化比特?；静襟E是：先用wavread 函數(shù)將WAV 文件轉(zhuǎn)換成列數(shù)組變量，再求出數(shù)組變量的極值并對所有元素作歸一化處理，最后用wavwrite 函數(shù)還原成音量標準化的WAV 文件。本文現(xiàn)對手機短信鈴聲的音頻信號進行處理。先將其復(fù)制另存到MATLAB 當前目錄中，文件名為test.wav 。再通過音量標準化處理后保存為outte

10、st.wav 文件。試聽可知標準化處理后音量稍大,在MATLAB運行結(jié)果中采樣頻率FS=44100, 量化比特NBITS=16，Ym=1.3.2聲道分離合并與組合立體聲或雙聲道音頻信號有左右兩個聲道。利用MATLAB 實現(xiàn)雙聲道分離， 兩路聲道合并和兩個單聲道組合成一個雙聲道等效果。實際上是利用了MATLAB 的矩陣抽取，矩陣相加和矩陣重組運算現(xiàn)以音量標準化處理生成的outtest.wav 實現(xiàn)分離，合并和組合處理?？梢栽嚶犅暤婪蛛x，合并與組合的效果。 也可對各文件大小進行比較。3.3數(shù)字濾波數(shù)字濾波是常用的音頻處理技術(shù)，可根據(jù)技術(shù)指標先利用FDAToo

11、l 工具設(shè)計一個數(shù)字濾波器，再用Filter 或Filter2 函數(shù)即可實現(xiàn)濾波處理。調(diào)用的Filter 函數(shù)格式是Y = filte (B,A,X) 。其中 B 和A 是濾波器傳輸函數(shù)的分子和分母系數(shù),X 是輸入變量，Y是實現(xiàn)濾波后的輸出變量。如果處理立體聲音頻信號，可分開處理，但用FIR 濾波器時調(diào)用Filter2 函數(shù)更方便。現(xiàn)以聲道分離合并與組合生成的outtest12.wav 實現(xiàn)數(shù)字濾波。程序運行結(jié)果如圖3所示，由圖可知濾波對波形影響不大，但對高頻有較大衰減，試聽

12、會感覺到處理后的聲音比較沉悶。圖3數(shù)字濾波前后的波形對比3.4數(shù)據(jù)轉(zhuǎn)換數(shù)據(jù)轉(zhuǎn)換是指改變音頻格式中的采樣頻率或量化位數(shù)。轉(zhuǎn)換原理是先用矩陣插值或抽取技術(shù)實現(xiàn)變量變換，如果是抽取數(shù)據(jù)還需在變換前作濾波處理使之滿足采樣定理。變量變換完成后再用Wavwrite 函數(shù)重新定義量化位數(shù)和采樣頻率即可實現(xiàn)數(shù)據(jù)轉(zhuǎn)換。數(shù)據(jù)轉(zhuǎn)換過程中要注意采樣頻率與原始采樣頻率及插值或抽取系數(shù)的關(guān)系。MATLAB 實現(xiàn)插值或抽取的函數(shù)有decimate，interp 和resample。這里以2 倍抽取為例，將經(jīng)過數(shù)字濾波后產(chǎn)生的outtest12Filter.WAV 文件進行數(shù)

13、據(jù)轉(zhuǎn)換處理。運行數(shù)據(jù)轉(zhuǎn)換程序，在得到的圖形窗口中執(zhí)行Edit/Axes Properties命令，再把各分圖下X 標簽中的Limits 設(shè)為0、0.01 和0、1000， 得到0-0.01 秒的波形和0 -1000Hz 的頻譜如圖4所示。由圖可知在滿足采樣定律條件下實現(xiàn)數(shù)據(jù)抽取。在原采樣率下，波形變密、頻譜變寬且幅度減半，但在新采樣率下波形和頻譜都很好。通過試聽輸出文件還可感受處理效果。圖4數(shù)據(jù)轉(zhuǎn)換波形對比圖音量標準化程序：clear; close all; clc;Y,FS,NBITS=wavread('test.WAV&#

14、39;); % 將WAV文件轉(zhuǎn)換成變量FS,NBITS, % 顯示采樣頻率和量化比特Ym=max(max(max(Y),max(abs(min(Y), % 找出雙聲道極值X=Y/Ym; % 歸一化處理wavwrite(X,FS,NBITS,'outtest.wav') % 將變量轉(zhuǎn)換成WAV文件聲道分離合并與組合程序：clear; close all; clc;x,FS,NBITS=wavread(' outtest.WAV'); % 將WAV文件轉(zhuǎn)換成變量x1=x(:,1); % 抽取第&#

15、160;1 聲道x2=x(:,2); % 抽取第 2 聲道wavwrite(x1,FS,NBITS,'outtest1.WAV'); % 實現(xiàn)1聲道分離wavwrite(x2,FS,NBITS,'outtest2.WAV'); % 實現(xiàn)2聲道分離%如果合并位置不對前面補0 %聲道長度不對后面補0x12=x1+x2; % 兩路單聲道列向量矩陣變量合并x12m=max(max(x12),abs(min(x12), % 找出極值y12=x12./x12m; % 歸一化處理wavw

16、rite(y12,FS,NBITS,'outtest12.WAV'); % 實現(xiàn)兩路聲道合并%如果組合位置不對前面補0-聲道長度不對后面補0x3=x1,x2; % 兩路單聲道變量組合wavwrite(x3,FS,NBITS,'outtest13.WAV'); % 實現(xiàn)兩路聲道組合數(shù)字濾波程序：clear; close all; clc;X,FS,NBITS=wavread('outtest12.WAV'); % 將WAV文件轉(zhuǎn)換成變量%利用 FDATool 設(shè)計一個 Lowpa

17、ssButterworth 濾波器%指標 FS=22050Hz Fp=1000Hz Ap=1dB Fs=3000Hz As=20dBB =0.0062,0.0187,0.0187,0.0062; % 分子系數(shù)A =1,-2.1706,1.6517,-0.4312; % 分母系數(shù)Y=filter(B,A,X); % 實現(xiàn)數(shù)字濾波t=(0:length(X)-1)/FS; % 計算數(shù)據(jù)時刻subplot(2,2,1);plot(t,X); % 繪制原波形圖title(' 原信號波形圖 ');

18、% 加標題subplot(2,2,3);plot(t,Y); % 繪制濾波波形圖title(' 濾波后波形圖 '); % 加標題xf=fft(X); % 作傅里葉變換求原頻譜yf=fft(Y); % 作傅里葉變換求濾波后頻譜fm=3000*length(xf)/FS; % 確定繪頻譜圖的上限頻率f=(0:fm)*FS/length(xf); % 確定繪頻譜圖的頻率刻度subplot(2,2,2);plot(f,abs(xf(1:length(f); % 繪制原波形頻譜圖title

19、(' 原信號頻譜圖 '); % 加標題subplot(2,2,4);plot(f,abs(yf(1:length(f); % 繪制濾波后頻譜圖title(' 濾波后信號頻譜圖 '); % 加標題wavwrite(Y,FS,NBITS,'outtest12Filter.WAV'); % 寫成WAV文件數(shù)據(jù)轉(zhuǎn)換程序：clear; close all; clc;x,FS,NBITS=wavread('outtestFilter.WAV'); % 將

20、WAV文件轉(zhuǎn)換成變量N=length(x); % 計算數(shù)據(jù)點數(shù)% 不是偶數(shù)點化成偶數(shù)點if mod(N,2)=0; N=N; else x(N)=; N=N-1; end;% 原信號波形頻譜分析tx=(0:N-1)/FS; % 計算原信號數(shù)據(jù)點時刻subplot(3,2,1);plot(tx,x); % 繪制原信號波形title(' 原信號波形圖 '); % 加標題xf=fft(x); % 求原信號頻譜fx=(0:N/2)*FS/N; % 確定頻譜圖頻率刻度subplot(3,2

21、,2);plot(fx,abs(xf(1:N/2+1); % 繪制原信號頻譜title(' 原信號頻譜圖 '); % 加標題% 實現(xiàn)數(shù)據(jù)抽取k=1:N/2; % 確定抽取位置y=x(2*k); % 實現(xiàn)抽取后的數(shù)據(jù)M=length(y); % 計算抽取后數(shù)據(jù)點數(shù)% 抽取數(shù)據(jù)在原采樣頻率 FS 下的波形頻譜分析ty=(0:M-1)/FS; % 計算數(shù)據(jù)點時刻subplot(3,2,3);plot(ty,y); % 繪制信號波形圖title(' 原采樣率下新波形圖 '); % 加標題yf=fft(y); % 求頻譜fy=(0:M/2)*FS/M; % 確定頻譜圖頻率刻度subplot(3,2,4);plot(fy,abs(yf(1:M/2+