基于MATLAB算法的數(shù)字變聲器.doc

基于MATLAB算法的數(shù)字變聲器【摘要】變聲器的原理是通過改變輸入聲音頻率，進行改變聲音的音色、音調(diào)，使輸出聲音在感官上與原聲音不同。變聲器是借助對聲音音色和音調(diào)的雙重復(fù)合改變，實現(xiàn)輸出聲音的改變。通過自己發(fā)聲，共振峰頻率的改變是基于重采樣實現(xiàn)的。同時用LPC倒譜分析法分析共振峰的范圍，通過MATLAB編寫程序并構(gòu)置GUI界面。實驗中通過MATLAB軟件對采集來的語言信號進行頻譜分析，討論不同人之間的噪音源參數(shù)和共振峰參數(shù)，通過對不同人語言信號頻率和幅度的改變來實現(xiàn)不同人之間的語音轉(zhuǎn)換?！娟P(guān)鍵詞】短時自相關(guān)法； LPC倒譜； 語音信號； matlab GUI；1引言隨著生活水平的提高，科技的不斷進步，很多人為了娛樂，從而希望改變自己的聲音；還有如今的許多的訪問節(jié)目為了保護被訪問者，都對聲音進行了相應(yīng)的處理。本設(shè)計通過編寫matlab程序，借助對聲音音色和音調(diào)的雙重復(fù)合改變，實現(xiàn)輸出聲音的改變。共振峰頻率的改變是基于重采樣實現(xiàn)的，從重采樣原理知道，這也同時印發(fā)了基頻的變化，為保證基頻變化和共振峰頻率變化的獨立、互不相關(guān)，在基頻移動時必須考慮抵消重采樣帶來的偏移，理論上只要基頻檢測足夠精確，確實可以保證基頻改變和共振峰頻率改變間的互不相關(guān)。保證變聲效果的自然度主要是沒有采用基因檢測將基因移動和共振峰變化徹底隔離的緣故。本次課程設(shè)計就是運用我們所學(xué)到的理論知識，用MATLAB軟件來實現(xiàn)對語音信號的變聲處理，理論聯(lián)系實際，從而更好地掌握以及運用所學(xué)習(xí)的知識。2數(shù)字變聲器的原理與算法2.1基本原理語音科學(xué)家將人類發(fā)聲過程視作一個由聲門源輸送的氣流經(jīng)以聲道、口、鼻腔組成的濾波器調(diào)制而成的。人類語音可分為有聲語音和無聲語音，前者是由聲帶振動激勵的脈沖信號經(jīng)聲腔調(diào)制變成不同的音，它是人類語言中元音的基礎(chǔ)，聲帶振動的頻率稱為基頻。無聲語音則是聲帶保持開啟狀態(tài)，禁止振動引發(fā)的。一般來說，由聲門振動決定的基頻跟說話人的性別特征有關(guān)，如下表，而無聲語音則沒有體現(xiàn)這個特征。說話人的個性化音色和語音的另外一個聲學(xué)參數(shù)共振峰頻率的分布有關(guān)。兒童由于聲道短，其共振峰頻率高于成年人，成年女性的聲道一般短于成年男性，所以女性的共振峰頻率一般高于男性。在進行性別變聲時，主要考慮基音周期、基頻和共振峰頻率的變化。其中男生、女生和童聲的基頻、共振峰的關(guān)系如圖1所示；基音周期改變時，基頻、共振峰同時變化，若伸展既有男變女、女變童，反之亦可。本實驗是基于打開一種聲音進行相關(guān)參數(shù)提取，修改接近于女聲、男聲或童聲，實現(xiàn)聲音的變換。 圖1 人群基頻與共振峰的關(guān)系2.2實現(xiàn)過程及算法采用線性預(yù)測參數(shù)合成法。線性預(yù)測參數(shù)合成法利用LPC語音分析方法，通過分析自然語音樣本，計算出LPC系數(shù)，就可以建立信號產(chǎn)生模型，從而合成出語音。大致流程如圖2所示。圖2 制作流程 2.2.1基于短時自相關(guān)法的基音周期估值對語言信號進行低通濾波，然后進行自相關(guān)計算。在低通濾波時，采用巴特沃斯濾波器。根據(jù)人的說話特征設(shè)定相應(yīng)指標(biāo)參數(shù)，對本段語音設(shè)計算出巴特沃斯模擬濾波器的階數(shù)N為5，3dB截止頻率c, c=(100.1ap-1)-12n,算出c為0.175，歸一化低通原型系統(tǒng)函數(shù)為：Ga（p）=1p5+b4p4+b3p3+b2p2+b1p+b0其中b0=1.0000,b1=3.2361,b2=5.2361,b3=5.2361,b4=3.2361將p=s/c帶人Ga（p）中，得到低通濾波器，Hs=c5s5+b4c2s4+b2c3s2+b1c4s+b0c5根據(jù)設(shè)定的濾波器編寫matlab程序，當(dāng)信號經(jīng)過低通濾波器后，對原始信號濾波產(chǎn)生結(jié)果如下圖2所示，低通濾波后，保留基音頻率，然后再用8kHz采樣頻率進行采樣，采樣序列為x(n)，然后進行下一步的自相關(guān)計算。圖3 原始信號低通濾波 2.2.2語音信號的短時自相關(guān)函數(shù)定義語音信號自相關(guān)函數(shù)如下：Rnk=n=0n=k-n-1xnnxn(n+k)其中k為信號延遲點數(shù): xnn為語音信號；N為語音幀長度。經(jīng)過低通濾波之后，取160個樣點數(shù)，幀長取10ms，對每幀語音求短時自相關(guān)，取得自相關(guān)最大點數(shù)，自相關(guān)函數(shù)在基音周期處表現(xiàn)為峰值，這些峰值點之間的間隔的平均值就是基音周期。 2.2.3 LPC倒譜法提取共振峰通過線性預(yù)測分析得到合成濾波器的系統(tǒng)函數(shù)為HZ=11-i=1paiz-i 其沖擊響應(yīng)為h(n),ai為預(yù)測系數(shù)。然后求h(n)的倒譜(n),首先根據(jù)同態(tài)分析方法，有(z)=logH(z),因為H（z）是最小相位的，即在單位圓內(nèi)是解析的，所以(z)可以展開成級數(shù)形式，即(z)=n=1(n)z-n,也就是說(z)的逆變換(n)是存在的，設(shè)(n)=0，將等式兩端分別對z-1求導(dǎo)，得到： n=1n(n)z-n+1=i=1piz-i+11-i=1paiz-i有（1-i=1paiz-i）n=1n(n)z-n+1=n=1iz-i+1, 令式左右兩邊Z的各次冪前系數(shù)分別相等得到：(1)=1(n)=an+i=1n-1(1-in)ai(n-i) 1np(n)=i=1n-1(1-in)ai(n-i) np在本實驗中，取預(yù)測階數(shù)為10，語音的樣點數(shù)為160，按照上式可直接從預(yù)測系數(shù)ai求得倒譜(n)，這個倒譜是根據(jù)線性預(yù)測模型得到，即稱為LPC倒譜。2.2.4線性預(yù)測語音信號合成根據(jù)線性預(yù)測的基本思想，用過去M個樣點值來預(yù)測現(xiàn)在或未來的樣點值：（n）=i=1Maiy(n-i)（n）是預(yù)測信號，ai為預(yù)測系數(shù)，y(n)為取樣信號，預(yù)測誤差（n）:E2(n)=Eyn-i=1Maiy(n-i)2為使E2(n)最小，對ai求偏導(dǎo)，并令其為零，有：Eyn-i=1Maiy(n-i)y(n-j)=0 j=1,M上式表明采用最佳預(yù)測系數(shù)時，預(yù)測誤差與過去的樣點值正交。由于語音信號的短時平穩(wěn)性，要分幀處理，對于每一個樣點值記為 yn(n), yn(n+N-1)，這段語音記為Yn。對于語音段Yn，并記n(j,i)=Eynm-jym(m-i),對于語音段Yn，它的自相關(guān)函數(shù)為： Rnj=m=jn-1ynmyn(m-j) n=1,M根據(jù)Yulewalker方程，可以解出樣值，用這種方法定期地改變激勵參數(shù)u(n)和預(yù)測系數(shù)ai，并使用修改過后的基因周期和共振峰參數(shù)，就能合成出語音，合成語音樣本如下：S（n）=i=1pais(n-i)+Gu(n)ai為預(yù)測系數(shù)；G為模型增益；u(n)為激勵；合成樣本s(n)；p為預(yù)測器階數(shù)；運行算法2.2.5流程圖開始開始打開音頻導(dǎo)入語音數(shù)據(jù)S 數(shù)據(jù)歸一化變童聲女變男男變女建立濾波器重建語言合成并加窗確定幀數(shù)線性預(yù)測法預(yù)測p個系數(shù)計算減小基音周期提取共振峰合成語音合成激勵3數(shù)字變聲器的軟件實現(xiàn)及仿真結(jié)果3.1 界面設(shè)計在設(shè)計數(shù)字變聲器的界面時，使用了MATLAB提供的可視化的界面環(huán)境Guide。相比通過編寫程序進行GUI的設(shè)計，可視化的界面環(huán)境更加方便快捷。見圖4圖4 界面設(shè)計 1按鈕組：設(shè)計中使用了兩個按鈕組，把打開、原聲、男變女、女變男、童聲三項功能放在一組；另一個把語音的圖像 即4個axes放在了一組。 2按鈕：最常用的控件，用于相應(yīng)用戶的鼠標(biāo)單擊，按鈕上有說明文字說明其作用。設(shè)計中使用了五個按鈕，分別用于打開、原聲、男變女、女變男、童聲。3 軸： 用于顯示圖像。3.2回調(diào)函數(shù)對控件的Callback屬性編程是實現(xiàn)GUI的基本機制，我們下來就要在各函數(shù)中編寫程序代碼，完成各回調(diào)函數(shù)的功能。 （1）打開function dakai_Callback(hObject, eventdata, handles)filename,pathname=uigetfile(*.wav,ALL FILES,選擇聲音);ifisequal(filename pathname,0,0)return;endstr=pathname filename;temp,Fs=audioread(str);temp1=resample(temp,80,441);handles.y=temp1;handles.y1=temp;handles.Fs=Fs;guidata(hObject,handles);（2）原聲function yuansheng_Callback(hObject, eventdata, handles)c=handles.Fs;sound(handles.y1,c);plot(handles.axes1,handles.y1)title(handles.axes1,時域圖);ysize=size(handles.y1);y1=fft(handles.y,length(handles.y1);ysize=size(y1);plot(handles.axes2,abs(y1);xlabel(handles.axes2,頻率);ylabel(handles.axes2,幅度);title(handles.axes2,頻率特性);（3）男變女function nanbiannv_Callback(hObject, eventdata, handles)FL = 80 ; % 幀移 WL = 240 ; % 窗長 P = 10 ; %預(yù)測系數(shù)個數(shù) s = handles.y;fs = handles.Fs;% 定義常數(shù) s = s/max(s); % 歸一化 L = length(s); % 讀入語音長度 FN = floor(L/FL)-2; % 計算幀長，floor；向負(fù)無窮方向% 預(yù)測和重建濾波器exc = zeros(L,1); % 激勵信號，double類零矩陣L行1列zi_pre = zeros(P,1); % 預(yù)測濾波器狀態(tài)s_rec = zeros(L,1); % 重建語音zi_rec = zeros(P,1);% 變調(diào)濾波器exc_syn_t = zeros(L,1); % 合成的激勵信號，創(chuàng)建一個L行1列的0脈沖s_syn_t = zeros(L,1); % 合成語音last_syn_t = 0; % 存儲上一個段的最后一個脈沖的下標(biāo)zi_syn_t = zeros(P,1); % 合成濾波器hw = hamming(WL); %漢明窗%濾波器% 依次處理每幀語音for n = 3:FN %從第三個子數(shù)組開始% 計算預(yù)測系數(shù)s_w = s(n*FL-WL+1:n*FL).*hw; %漢明窗加權(quán) A,E=lpc(s_w,P); %線性預(yù)測計算預(yù)測系數(shù)% A是預(yù)測系數(shù)，E會被用來計算合成激勵的能量s_f=s(n-1)*FL+1:n*FL); % 本幀語音%利用filter函數(shù)重建語音 exc1,zi_pre = filter(A,1,s_f,zi_pre); exc(n-1)*FL+1:n*FL) = exc1; %計算激勵%利用filter函數(shù)重建語音 s_rec1,zi_rec = filter(1,A,exc1,zi_rec);s_rec(n-1)*FL+1:n*FL) = s_rec1; %重建語音% 下面只有得到exc后才可以s_Pitch = exc(n*FL-222:n*FL); PT(n) = findpitch(s_Pitch); %計算基音周期pt G = sqrt(E*PT(n); %計算合成激勵的能量G% tempn_syn = 1:n*FL-last_syn;% exc_syn1 = zeros(length(tempn_syn),1);% exc_syn1(mod(tempn_syn,PT)=0) = G; %某一段算出的脈沖% exc_syn1 = exc_syn1(n-1)*FL-last_syn+1:n*FL-last_syn);% s_syn1,zi_syn = filter(1,A,exc_syn1,zi_syn);% exc_syn(n-1)*FL+1:n*FL)=exc_syn1; %計算得到合成激勵% s_syn(n-1)*FL+1:n*FL) = s_syn1; %計算得到合成語音% last_syn = last_syn+PT*floor(n*FL-last_syn)/PT); PT1 =floor(PT(n)/2); %減小基音周期poles = roots(A);deltaOMG =100*2*pi/fs;for p=1:10 %增加共振峰ifimag(poles(p)0poles(p) = poles(p)*exp(1j*deltaOMG);elseifimag(poles(p)=0.85*Rop Rop=R2; Top=T2;endif R30.85*Rop Rop=R3; Top=T3;endPT=Top;return（3）女變男 function nvbiannan_Callback(hObject, eventdata, handles)FL = 80 ; % 幀移 WL = 240 ; % 窗長 P = 10 ; %預(yù)測系數(shù)個數(shù) s = handles.y;fs = handles.Fs;% 定義常數(shù) s = s/max(s); % 歸一化 L = length(s); % 讀入語音長度 FN = floor(L/FL)-2; % 計算幀長，floor；向負(fù)無窮方向% 預(yù)測和重建濾波器exc = zeros(L,1); % 激勵信號，double類零矩陣L行1列zi_pre = zeros(P,1); % 預(yù)測濾波器狀態(tài)s_rec = zeros(L,1); % 重建語音zi_rec = zeros(P,1);% 變調(diào)濾波器exc_syn_t = zeros(L,1); % 合成的激勵信號，創(chuàng)建一個L行1列的0脈沖s_syn_t = zeros(L,1); % 合成語音last_syn_t = 0; % 存儲上一個段的最后一個脈沖的下標(biāo)zi_syn_t = zeros(P,1); % 合成濾波器hw = hamming(WL); %漢明窗%濾波器% 依次處理每幀語音for n = 3:FN %從第三個子數(shù)組開始% 計算預(yù)測系數(shù)s_w = s(n*FL-WL+1:n*FL).*hw; %漢明窗加權(quán) A,E=lpc(s_w,P); %線性預(yù)測計算預(yù)測系數(shù)% A是預(yù)測系數(shù)，E會被用來計算合成激勵的能量s_f=s(n-1)*FL+1:n*FL); % 本幀語音%利用filter函數(shù)重建語音 exc1,zi_pre = filter(A,1,s_f,zi_pre); exc(n-1)*FL+1:n*FL) = exc1; %計算激勵%利用filter函數(shù)重建語音 s_rec1,zi_rec = filter(1,A,exc1,zi_rec);s_rec(n-1)*FL+1:n*FL) = s_rec1; %重建語音% 下面只有得到exc后才可以s_Pitch = exc(n*FL-222:n*FL); PT(n) = findpitch(s_Pitch); %計算基音周期pt G = sqrt(E*PT(n); %計算合成激勵的能量G% tempn_syn = 1:n*FL-last_syn;% exc_syn1 = zeros(length(tempn_syn),1);% exc_syn1(mod(tempn_syn,PT)=0) = G; %某一段算出的脈沖% exc_syn1 = exc_syn1(n-1)*FL-last_syn+1:n*FL-last_syn);% s_syn1,zi_syn = filter(1,A,exc_syn1,zi_syn);% exc_syn(n-1)*FL+1:n*FL)=exc_syn1; %計算得到合成激勵% s_syn(n-1)*FL+1:n*FL) = s_syn1; %計算得到合成語音% last_syn = last_syn+PT*floor(n*FL-last_syn)/PT); PT1 =floor(PT(n)/2); %減小基音周期poles = roots(A);deltaOMG =100*2*pi/fs;for p=1:10 %增加共振峰ifimag(poles(p)0poles(p) = poles(p)*exp(1j*deltaOMG);elseifimag(poles(p)0 ; poles(p) = poles(p)*exp(1i*deltaOMG); else if imag(poles(p)0 ; poles(p) = poles(p)*exp(-1i*deltaOMG); end endend A1=poly(poles); tempn_syn_t = 1:n*FL-last_syn_t; exc_syn1_t = zeros(length(tempn_syn_t),1); exc_syn1_t(mod(tempn_syn_t,PT1)=0) = G; exc_syn1_t = exc_syn1_t(n-1)*FL-last_syn_t+1:n*FL-last_syn_t); s_syn1_t,zi_syn_t = filter(1,A1,exc_syn1_t,zi_syn_t); exc_syn_t(n-1)*FL+1:n*FL) = exc_syn1_t; s_syn_t(n-1)*FL+1:n*FL) = s_syn1_t; last_syn_t = last_syn_t+PT1*floor(n*FL-last_syn_t)/PT1);end plot(handles.axes4,exc_syn_t) xlabel(handles.axes4,時域圖); ylabel(handles.axes4,頻率); title(handles.