基于MATLAB算法的数字变声器

1、基于MATLAB算法的数字变声器【摘要】变声器的原理是通过改变输入声音频率，进行改变声音的音色、音调，使输出声音在感官上与原声音不同。变声器是借助对声音音色和音调的双重复合改变，实现输出声音的改变。通过自己发声，共振峰频率的改变是基于重采样实现的。同时用LPC倒谱分析法分析共振峰的范围，通过MATLAB编写程序并构置GUI界面。实验中通过MATLAB软件对采集来的语言信号进行频谱分析，讨论不同人之间的噪音源参数和共振峰参数，通过对不同人语言信号频率和幅度的改变来实现不同人之间的语音转换。【关键词】短时自相关法； LPC倒谱； 语音信号； matlab GUI；1引言随着生活水平的提高，科技的不

2、断进步，很多人为了娱乐，从而希望改变自己的声音；还有如今的许多的访问节目为了保护被访问者，都对声音进行了相应的处理。本设计通过编写matlab程序，借助对声音音色和音调的双重复合改变，实现输出声音的改变。共振峰频率的改变是基于重采样实现的，从重采样原理知道，这也同时印发了基频的变化，为保证基频变化和共振峰频率变化的独立、互不相关，在基频移动时必须考虑抵消重采样带来的偏移，理论上只要基频检测足够精确，确实可以保证基频改变和共振峰频率改变间的互不相关。保证变声效果的自然度主要是没有采用基因检测将基因移动和共振峰变化彻底隔离的缘故。本次课程设计就是运用我们所学到的理论知识，用MATLAB软件来实现对

3、语音信号的变声处理，理论联系实际，从而更好地掌握以及运用所学习的知识。2数字变声器的原理与算法2.1基本原理语音科学家将人类发声过程视作一个由声门源输送的气流经以声道、口、鼻腔组成的滤波器调制而成的。人类语音可分为有声语音和无声语音，前者是由声带振动激励的脉冲信号经声腔调制变成不同的音，它是人类语言中元音的基础，声带振动的频率称为基频。无声语音则是声带保持开启状态，禁止振动引发的。一般来说，由声门振动决定的基频跟说话人的性别特征有关，如下表，而无声语音则没有体现这个特征。说话人的个性化音色和语音的另外一个声学参数共振峰频率的分布有关。儿童由于声道短，其共振峰频率高于成年人，成年女性的声道一般短

4、于成年男性，所以女性的共振峰频率一般高于男性。在进行性别变声时，主要考虑基音周期、基频和共振峰频率的变化。其中男生、女生和童声的基频、共振峰的关系如图1所示；基音周期改变时，基频、共振峰同时变化，若伸展既有男变女、女变童，反之亦可。本实验是基于打开一种声音进行相关参数提取，修改接近于女声、男声或童声，实现声音的变换。 图1 人群基频与共振峰的关系2.2实现过程及算法采用线性预测参数合成法。线性预测参数合成法利用LPC语音分析方法，通过分析自然语音样本，计算出LPC系数，就可以建立信号产生模型，从而合成出语音。大致流程如图2所示。图2 制作流程 2.2.1基于短时自相关法的基音周期估值对语言信号

5、进行低通滤波，然后进行自相关计算。在低通滤波时，采用巴特沃斯滤波器。根据人的说话特征设定相应指标参数，对本段语音设计算出巴特沃斯模拟滤波器的阶数N为5，3dB截止频率c, c=(100.1ap-1)-12n,算出c为0.175，归一化低通原型系统函数为：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根据设定的滤波器编写matlab程序，当信号经过低通滤波器后，对原始信号

6、滤波产生结果如下图2所示，低通滤波后，保留基音频率，然后再用8kHz采样频率进行采样，采样序列为x(n)，然后进行下一步的自相关计算。图3 原始信号低通滤波 2.2.2语音信号的短时自相关函数定义语音信号自相关函数如下：Rnk=n=0n=k-n-1xnnxn(n+k)其中k为信号延迟点数: xnn为语音信号；N为语音帧长度。经过低通滤波之后，取160个样点数，帧长取10ms，对每帧语音求短时自相关，取得自相关最大点数，自相关函数在基音周期处表现为峰值，这些峰值点之间的间隔的平均值就是基音周期。 2.2.3 LPC倒谱法提取共振峰通过线性预测分析得到合成滤波器的系统函数为HZ=11-i=1pai

7、z-i 其冲击响应为h(n),ai为预测系数。然后求h(n)的倒谱(n),首先根据同态分析方法，有(z)=logH(z),因为H（z）是最小相位的，即在单位圆内是解析的，所以(z)可以展开成级数形式，即(z)=n=1(n)z-n,也就是说(z)的逆变换(n)是存在的，设(n)=0，将等式两端分别对z-1求导，得到： 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的各次幂前系数分别相等得到：(1)=1(n)=an+i=1n-1(1-in)ai(n-i) 1np(n)=i=1n-1(1-

8、in)ai(n-i) n>p在本实验中，取预测阶数为10，语音的样点数为160，按照上式可直接从预测系数ai求得倒谱(n)，这个倒谱是根据线性预测模型得到，即称为LPC倒谱。2.2.4线性预测语音信号合成根据线性预测的基本思想，用过去M个样点值来预测现在或未来的样点值：（n）=i=1Maiy(n-i)（n）是预测信号，ai为预测系数，y(n)为取样信号，预测误差（n）:E2(n)=Eyn-i=1Maiy(n-i)2为使E2(n)最小，对ai求偏导，并令其为零，有：Eyn-i=1Maiy(n-i)y(n-j)=0 j=1,M上式表明采用最佳预测系数时，预测误差与过去的样点值正交。由于语音信

9、号的短时平稳性，要分帧处理，对于每一个样点值记为 yn(n), yn(n+N-1)，这段语音记为Yn。对于语音段Yn，并记n(j,i)=Eynm-jym(m-i),对于语音段Yn，它的自相关函数为： Rnj=m=jn-1ynmyn(m-j) n=1,M根据Yulewalker方程，可以解出样值，用这种方法定期地改变激励参数u(n)和预测系数ai，并使用修改过后的基因周期和共振峰参数，就能合成出语音，合成语音样本如下：S（n）=i=1pais(n-i)+Gu(n)ai为预测系数；G为模型增益；u(n)为激励；合成样本s(n)；p为预测器阶数；运行算法2.2.5流程图开始开始打开音频导入语音数据S

10、 数据归一化变童声女变男男变女建立滤波器重建语言合成并加窗确定帧数线性预测法预测p个系数计算减小基音周期提取共振峰合成语音合成激励3数字变声器的软件实现及仿真结果3.1 界面设计在设计数字变声器的界面时，使用了MATLAB提供的可视化的界面环境Guide。相比通过编写程序进行GUI的设计，可视化的界面环境更加方便快捷。见图4图4 界面设计 1按钮组：设计中使用了两个按钮组，把打开、原声、男变女、女变男、童声三项功能放在一组；另一个把语音的图像 即4个axes放在了一组。 2按钮：最常用的控件，用于相应用户的鼠标单击，按钮上有说明文字说明其作用。设计中使用了五个按钮，分别用于打开、原声、男变女、

11、女变男、童声。3 轴： 用于显示图像。3.2回调函数对控件的Callback属性编程是实现GUI的基本机制，我们下来就要在各函数中编写程序代码，完成各回调函数的功能。 （1）打开function dakai_Callback(hObject, eventdata, handles)filename,pathname=uigetfile('*.wav','ALL FILES','选择声音');ifisequal(filename pathname,0,0)return;endstr=pathname filename;temp,Fs=audiore

12、ad(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(hand

13、les.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 ; %预测系数个数 s = handles.y;fs

14、= handles.Fs;% 定义常数 s = s/max(s); % 归一化 L = length(s); % 读入语音长度 FN = floor(L/FL)-2; % 计算帧长，floor；向负无穷方向% 预测和重建滤波器exc = zeros(L,1); % 激励信号，double类零矩阵L行1列zi_pre = zeros(P,1); % 预测滤波器状态s_rec = zeros(L,1); % 重建语音zi_rec = zeros(P,1);% 变调滤波器exc_syn_t = zeros(L,1); % 合成的激励信号，创建一个L行1列的0脉冲s_syn_t = zeros(L,1

15、); % 合成语音last_syn_t = 0; % 存储上一个段的最后一个脉冲的下标zi_syn_t = zeros(P,1); % 合成滤波器hw = hamming(WL); %汉明窗%滤波器% 依次处理每帧语音for n = 3:FN %从第三个子数组开始% 计算预测系数s_w = s(n*FL-WL+1:n*FL).*hw; %汉明窗加权 A,E=lpc(s_w,P); %线性预测计算预测系数% A是预测系数，E会被用来计算合成激励的能量s_f=s(n-1)*FL+1:n*FL); % 本帧语音%利用filter函数重建语音 exc1,zi_pre = filter(A,1,s_f,

16、zi_pre); exc(n-1)*FL+1:n*FL) = exc1; %计算激励%利用filter函数重建语音 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(leng

17、th(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)/P

18、T); 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(-1j*deltaOMG);endendend A1=poly(poles);tempn_syn_t=(1:n*FL-last_syn_t); exc_syn1_t = zeros(length(tempn_sy

19、n_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

20、)/PT1);end%绘图plot(handles.axes4,exc_syn_t) xlabel(handles.axes4,'时域图'); ylabel(handles.axes4,'频率'); title(handles.axes4,'幅度'); handles.y=s_syn_t; guidata(hObject,handles); plot(handles.axes3,s_syn_t); %y1size=size(s_syn_t); plot(handles.axes3,s_syn_t); title(handles.axes3,

21、9;频率特性'); sound(handles.y,8000);%findpitch函数 function PT=findpitch(s)B,A=butter(5,700/4000);s=filter(B,A,s);R=zeros(143,1);for k=1:143 R(k)=s(144:223)'*s(144-k:223-k);endR1,T1=max(R(80:143);T1=T1+79;R1=R1/(norm(s(144-T1:223-T1)+1);R2,T2=max(R(40:79);T2=T2+39;R2=R2/(norm(s(144-T2:223-T2)+1);R

22、3,T3=max(R(20:39);T3=T3+19;R3=R3/(norm(s(144-T3:223-T3)+1);Top=T1;Rop=R1;if R2>=0.85*Rop Rop=R2; Top=T2;endif R3>0.85*Rop Rop=R3; Top=T3;endPT=Top;return（3）女变男 function nvbiannan_Callback(hObject, eventdata, handles)FL = 80 ; % 帧移 WL = 240 ; % 窗长 P = 10 ; %预测系数个数 s = handles.y;fs = handles.Fs;

23、% 定义常数 s = s/max(s); % 归一化 L = length(s); % 读入语音长度 FN = floor(L/FL)-2; % 计算帧长，floor；向负无穷方向% 预测和重建滤波器exc = zeros(L,1); % 激励信号，double类零矩阵L行1列zi_pre = zeros(P,1); % 预测滤波器状态s_rec = zeros(L,1); % 重建语音zi_rec = zeros(P,1);% 变调滤波器exc_syn_t = zeros(L,1); % 合成的激励信号，创建一个L行1列的0脉冲s_syn_t = zeros(L,1); % 合成语音last

24、_syn_t = 0; % 存储上一个段的最后一个脉冲的下标zi_syn_t = zeros(P,1); % 合成滤波器hw = hamming(WL); %汉明窗%滤波器% 依次处理每帧语音for n = 3:FN %从第三个子数组开始% 计算预测系数s_w = s(n*FL-WL+1:n*FL).*hw; %汉明窗加权 A,E=lpc(s_w,P); %线性预测计算预测系数% A是预测系数，E会被用来计算合成激励的能量s_f=s(n-1)*FL+1:n*FL); % 本帧语音%利用filter函数重建语音 exc1,zi_pre = filter(A,1,s_f,zi_pre); exc(

25、n-1)*FL+1:n*FL) = exc1; %计算激励%利用filter函数重建语音 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)

26、,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 =floo

27、r(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(-1j*deltaOMG);endendend A1=poly(poles);tempn_syn_t=(1:n*FL-last_syn_t); exc_syn1_t = zeros(length(tempn_syn_t),1); exc_

28、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%绘图

29、handles.y=s_syn_t;guidata(hObject,handles);%绘图plot(handles.axes4,exc_syn_t) xlabel(handles.axes4,'时域图'); ylabel(handles.axes4,'频率'); title(handles.axes4,'幅度'); handles.y=s_syn_t; guidata(hObject,handles); plot(handles.axes3,s_syn_t); %y1size=size(s_syn_t); plot(handles.axes3,

30、s_syn_t); title(handles.axes3,'频率特性');（4）变童声function tongsheng_Callback(hObject, eventdata, handles)FL = 80 ; % 帧移 WL = 240 ; % 窗长 P = 10 ; %预测系数个数 s = handles.y;fs = handles.Fs;% 定义常数 s = s/max(s); % 归一化 L = length(s); % 读入语音长度 FN = floor(L/FL)-2; % 计算帧长，floor；向负无穷方向% 预测和重建滤波器exc = zeros(L,

31、1); % 激励信号，double类零矩阵L行1列zi_pre = zeros(P,1); % 预测滤波器状态s_rec = zeros(L,1); % 重建语音zi_rec = zeros(P,1);% 变调滤波器exc_syn_t = zeros(L,1); % 合成的激励信号，创建一个L行1列的0脉冲s_syn_t = zeros(L,1); % 合成语音last_syn_t = 0; % 存储上一个段的最后一个脉冲的下标zi_syn_t = zeros(P,1); % 合成滤波器hw = hamming(WL); %汉明窗%滤波器% 依次处理每帧语音for n = 3:FN %从第三个

32、子数组开始% 计算预测系数s_w = s(n*FL-WL+1:n*FL).*hw; %汉明窗加权 A,E=lpc(s_w,P); %线性预测计算预测系数% A是预测系数，E会被用来计算合成激励的能量s_f=s(n-1)*FL+1:n*FL); % 本帧语音%利用filter函数重建语音 exc1,zi_pre = filter(A,1,s_f,zi_pre); exc(n-1)*FL+1:n*FL) = exc1; %计算激励%利用filter函数重建语音 s_rec1,zi_rec = filter(1,A,exc1,zi_rec);s_rec(n-1)*FL+1:n*FL) = s_rec

33、1; %重建语音% 下面只有得到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

34、,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/2); poles = roots(A); deltaOMG = 700*2*pi/8000; for p=1:10 if imag(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+