基于MATLAB算法数字变声器

基于MATLAB算法的数字变声器【摘要】变声器的原理是通过改变输入声音频率，进行改变声音的音色、音调，使输出声音在感官上与原声音不同。

变声器是借助对声音音色和音调的双重复合改变，实现输出声音的改变。

通过自己发声，共振峰频率的改变是基于重采样实现的。

同时用LPC倒谱分析法分析共振峰的范围，通过MATLAB编写程序并构置GUI界面。

实验中通过MATLAB软件对采集来的语言信号进行频谱分析，讨论不同人之间的噪音源参数和共振峰参数，通过对不同人语言信号频率和幅度的改变来实现不同人之间的语音转换。

【关键词】短时自相关法；LPC倒谱；语音信号；matlab GUI；1引言随着生活水平的提高，科技的不断进步，很多人为了娱乐，从而希望改变自己的声音；还有如今的许多的访问节目为了保护被访问者，都对声音进行了相应的处理。

本设计通过编写matlab程序，借助对声音音色和音调的双重复合改变，实现输出声音的改变。

共振峰频率的改变是基于重采样实现的，从重采样原理知道，这也同时印发了基频的变化，为保证基频变化和共振峰频率变化的独立、互不相关，在基频移动时必须考虑抵消重采样带来的偏移，理论上只要基频检测足够精确，确实可以保证基频改变和共振峰频率改变间的互不相关。

保证变声效果的自然度主要是没有采用基因检测将基因移动和共振峰变化彻底隔离的缘故。

本次课程设计就是运用我们所学到的理论知识，用MATLAB软件来实现对语音信号的变声处理，理论联系实际，从而更好地掌握以及运用所学习的知识。

2数字变声器的原理与算法2.1基本原理语音科学家将人类发声过程视作一个由声门源输送的气流经以声道、口、鼻腔组成的滤波器调制而成的。

人类语音可分为有声语音和无声语音，前者是由声带振动激励的脉冲信号经声腔调制变成不同的音，它是人类语言中元音的基础，声带振动的频率称为基频。

无声语音则是声带保持开启状态，禁止振动引发的。

一般来说，由声门振动决定的基频跟说话人的性别特征有关，如下表，而无声语音则没有体现这个特征。

说话人的个性化音色和语音的另外一个声学参数——共振峰频率的分布有关。

儿童由于声道短，其共振峰频率高于成年人，成年女性的声道一般短于成年男性，所以女性的共振峰频率一般高于男性。

在进行性别变声时，主要考虑基音周期、基频和共振峰频率的变化。

其中男生、女生和童声的基频、共振峰的关系如图1所示；基音周期改变时，基频、共振峰同时变化，若伸展既有男变女、女变童，反之亦可。

本实验是基于打开一种声音进行相关参数提取，修改接近于女声、男声或童声，实现声音的变换。

2.2实现过程及算法采用线性预测参数合成法。

线性预测参数合成法利用LPC语音分析方法，通过分析自然语音样本，计算出LPC系数，就可以建立信号产生模型，从而合成出语音。

大致流程如图2所示。

2.2.1基于短时自相关法的基音周期估值对语言信号进行低通滤波，然后进行自相关计算。

在低通滤波时，采用巴特沃斯滤波器。

根据人的说话特征设定相应指标参数，对本段语音设计算出巴特沃斯模拟滤波器的阶数N为5，3dB截止频率Ωc,Ωc=(100.1a p−1)−12n,算出Ωc为0.175，归一化低通原型系统函数为：G a（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带人G a（p）中，得到低通滤波器，H s=Ωc5s5+b4Ωc2s4+b2Ωc3s2+b1Ωc4s+b0Ωc5根据设定的滤波器编写matlab程序，当信号经过低通滤波器后，对原始信号滤波产生结果如下图2所示，低通滤波后，保留基音频率，然后再用8kHz采样频率进行采样，采样序列为x(n)，然后进行下一步的自相关计算。

2.2.2语音信号的短时自相关函数图3 原始信号低通滤波图2 制作流程定义语音信号自相关函数如下：R n (k )=∑x n (n )x n (n +k )n =k −n −1n =0其中k 为信号延迟点数: x n (n )为语音信号；N 为语音帧长度。

经过低通滤波之后，取160个样点数，帧长取10ms ，对每帧语音求短时自相关，取得自相关最大点数，自相关函数在基音周期处表现为峰值，这些峰值点之间的间隔的平均值就是基音周期。

2.2.3 LPC 倒谱法提取共振峰通过线性预测分析得到合成滤波器的系统函数为H (Z )=11−∑a i z −ipi =1 其冲击响应为h(n),a i 为预测系数。

然后求h(n)的倒谱h ˆ(n),首先根据同态分析方法，有Hˆ(z)=logH(z),因为H （z ）是最小相位的，即在单位圆内是解析的，所以H ˆ(z)可以展开成级数形式，即(z)=∑h ˆ(n )z −n ∞n =1,也就是说(z)的逆变换(n)是存在的，设(n)=0，将等式两端分别对z −1求导，得到：∑n hˆ(n )∞n =1z−n +1=∑ui z −i +1pi =11−∑a i z−ipi =1有（1-∑a i z −i pi =1）∑n (n )∝n =1z−n +1=∑u i z −i +1∝n =1, 令式左右两边Z 的各次幂前系数分别相等得到：(1)=1 (n)=a n +∑(1−i n )a i (n −i )n −1i =1 1≤n ≤p (n)=∑(1−i n)a i h ˆ(n −i )n −1i =1 n>p 在本实验中，取预测阶数为10，语音的样点数为160，按照上式可直接从预测系数a i 求得倒谱(n)，这个倒谱是根据线性预测模型得到，即称为LPC 倒谱。

2.2.4线性预测语音信号合成根据线性预测的基本思想，用过去M 个样点值来预测现在或未来的样点值：y ˆ（n ）=∑a i y (n −i )M i =1H ˆH ˆh ˆh ˆh ˆhˆh ˆh ˆh ˆhˆyˆ（n ）是预测信号，a i 为预测系数，y(n)为取样信号，预测误差（n ）: E ｛ε2(n)｝=E ｛[y (n )−∑a i y (n −i )M i =1]2｝为使E ｛ε2(n)｝最小，对a i 求偏导，并令其为零，有：E ｛［y (n )−∑a i y (n −i )M i =1］y(n-)｝=0 =1,…,M上式表明采用最佳预测系数时，预测误差与过去的样点值正交。

由于语音信号的短时平稳性，要分帧处理，对于每一个样点值记为 y n (n)…, y n (n+N-1)，这段语音记为Yn 。

对于语音段Yn ，并记Φn (j,i)=E ｛y n (m −j )y m (m −i )｝,对于语音段Yn ，它的自相关函数为：R n (j )=∑y n (m )y n (m −j )n −1m =j n=1,…,M 根据Yule —walker 方程，可以解出样值，用这种方法定期地改变激励参数u(n)和预测系数a i ，并使用修改过后的基因周期和共振峰参数，就能合成出语音，合成语音样本如下：S （n ）=∑a i s (n −i )pi =1+Gu(n)a i 为预测系数；G 为模型增益；u(n)为激励；合成样本s(n)；p 为预测器阶数；2.2.5流程图开始 导入语音数据S 算法 数据归一化 建立滤波器重建语言 合成并加窗 确定帧数线性预测法预测p 个系数 计算 减小基音周期 打开音频 男变女 开始运行 女变男 变童声3数字变声器的软件实现及仿真结果 3.1 界面设计在设计数字变声器的界面时，使用了MATLAB 提供的可视化的界面环境Guide 。

相比通 过编写程序进行GUI 的设计，可视化的界面环境更加方便快捷。

见图41按钮组：设计中使用了两个按钮组，把打开、原声、男变女、女变男、童声三项功能放在一组；另一个把语音的图像 即4个axes 放在了一组。

2按钮：最常用的控件，用于相应用户的鼠标单击，按钮上有说明文字说明其作用。

设计中使用了五个按钮，分别用于打开、原声、男变女、女变男、童声。

3 轴： 用于显示图像。

3.2回调函数对控件的Callback 属性编程是实现GUI 的基本机制，我们下来就要在各函数中编写程序代码，完成各回调函数的功能。

（1）打开function dakai_Callback(hObject, eventdata, handles)[filename,pathname]=uigetfile({'*.wav','ALL FILES'},'选择声音'); if isequal([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 ; %预测系数个数s = handles.y;fs = 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); % 合成语音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,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); %计算基音周期ptG = 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 %增加共振峰if imag(poles(p))>0poles(p) = poles(p)*exp(1j*deltaOMG);elseif imag(poles(p))<0poles(p) = poles(p)*exp(-1j*deltaOMG);endendendA1=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.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,'频率特性');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:143R(k)=s(144:223)'*s(144-k:223-k);end[R1,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);[R3,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*RopRop=R2;Top=T2;endif R3>0.85*RopRop=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;% 定义常数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_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((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); %计算基音周期ptG = 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 %增加共振峰if imag(poles(p))>0poles(p) = poles(p)*exp(1j*deltaOMG);elseif imag(poles(p))<0poles(p) = poles(p)*exp(-1j*deltaOMG);endendendA1=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%绘图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,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,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 %从第三个子数组开始% 计算预测系数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_rec1; %重建语音% 下面只有得到exc后才可以s_Pitch = exc(n*FL-222:n*FL);PT(n) = findpitch(s_Pitch); %计算基音周期ptG = 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/2);poles = roots(A);deltaOMG = 700*2*pi/8000;for p=1:10if imag(poles(p))>0 ;poles(p) = poles(p)*exp(1i*deltaOMG);elseif imag(poles(p))<0 ;poles(p) = poles(p)*exp(-1i*deltaOMG);endendendA1=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);endplot(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,'频率特性');sound(handles.y,8000)3.2.1仿真结果（1）男声变女声前后图形对（2）女生变男生前后图形对比（3）男声变童声图形对比图6男声变童声前后图形对比4 结束语在开始制作时毫无头绪，但经过老师和同学的讲解，及自己所找到的各种文献及例子一点一点的探索，最终一点点做出来了。