基于MATLAB算法的数字变声器

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

【关键词】短时自相关法； LPC倒谱；语音信号； matlab GUI；

1引言

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

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

2数字变声器的原理与算法

2.1基本原理

语音科学家将人类发声过程视作一个由声门源输送的气流经以声道、口、鼻腔组成的滤波器调制而成的。人类语音可分为有声语音和无声语音，前者是由声带振动激励的脉冲信号经声腔调制变成不同的音，它是人类语言中元音的基础，声带振动的频率称为基频。无声语音则是声带保持开启状态，禁止振动引发的。一般来说，由声门振动决定的基频跟说话人的性别特征有关，如下表，而无声语音则没有体现这个特征。说话人的个性化音色和语音的另外一个声学参数——共振峰频率的分布有关。儿童由于声道短，其共振峰频率高于成年人，成年女性的声道一般短于成年男性，所以女性的共振峰频率一般高于男性。

在进行性别变声时，主要考虑基音周期、基频和共振峰频率的变化。其中男生、女生和童声的基频、共振峰的关系如图1所示；基音周期改变时，基频、共振峰同时变化，若伸展既有男变女、女变童，反之亦可。本实验是基于打开一种声音进行相关参数提取，修改

接近于女声、男声或童声，实现声音的变换。

图1 人群基频与共振峰的关系

2.2实现过程及算法

采用线性预测参数合成法。线性预测参数合成法利用LPC语音分析方法，通过分析自然

语音样本，计算出LPC系数，就可以建立信号产生模型，从而合成出语音。大致流程如图2所示。

2.2.1基于短时自相关法的基音周期估值

对语言信号进行低通滤波，然后进行自相关计算。在低通滤波时，采用巴特沃斯滤波器。

根据人的说话特征设定相应指标参数，对本段语音设计算出巴特沃斯模拟滤波器的阶数N为5，3dB截止频率,,算出为0.175，归一化低通原型系统函数为：

其中

将带人中，得到低通滤波器，

根据设定的滤波器编写matlab程序，当信号经过低通滤波器后，对原始信号滤波产生结果如下图2所示，低通滤波后，保留基音频率，然后再用8kHz采样频率进行采样，采样序列为x(n)，然后进行下一步的自相关计算。

2.2.2语音信号的短时自相关函数

定义语音信号自相关函数如下：

图3 原始信号低通滤波

图2 制作流程

其中k 为信号延迟点数:

为语音信号；N 为语音帧长度。经过低通滤波之后，

取160个样点数，帧长取10ms ，对每帧语音求短时自相关，取得自相关最大点

数，自相关函数在基音周期处表现为峰值，这些峰值点之间的间隔的平均值就是基音周期。

2.2.3 LPC 倒谱法提取共振峰

通过线性预测分析得到合成滤波器的系统函数为

其冲击响应为h(n),为预测系数。然后求h(n)的倒谱h ?

(n),首先根据同态分析方

法，有H

?(z)=logH(z),因为H （z ）是最小相位的，即在单位圆内是解析的，所以H ?(z)可以展开成级数形式，即H ?(z)=,也就是说H

?(z)的逆变换h ?

(n)是存在的，设h

?(n)=0，将等式两端分别对求导，得到： =

有（1-）

=

, 令式左右两边Z 的各次幂前系数分别

相等得到：

h ?(1)=1 h ?(n)=+ 1

h

?(n)= n>p

在本实验中，取预测阶数为10，语音的样点数为160，按照上式可直接从预测系数求

得倒谱h

?(n)，这个倒谱是根据线性预测模型得到，即称为LPC 倒谱。 2.2.4线性预测语音信号合成

根据线性预测的基本思想，用过去M 个样点值来预测现在或未来的样点值：

y

?（n ）=

y

?（n ）是预测信号，为预测系数，y(n)为取样信号，预测误差（n ）: E ｛(n)｝=E ｛

｝

为使E｛(n)｝最小，对求偏导，并令其为零，有：

E｛［］y(n-)｝=0 =1,…,M

上式表明采用最佳预测系数时，预测误差与过去的样点值正交。由于语音信号的短时平稳性，要分帧处理，对于每一个样点值记为(n)…, (n+N-1)，这段语音记为Yn。

对于语音段Yn，并记(j,i)=E｛｝,对于语音段Yn，它的自相关函数为：

n=1,…,M

根据Yule—walker方程，可以解出样值，用这种方法定期地改变激励参数u(n)和预测系数，并使用修改过后的基因周期和共振峰参数，就能合成出语音，合成语音样本如下：S（n）=+Gu(n)

为预测系数；G为模型增益；u(n)为激励；合成样本s(n)；p为预测器阶数；

2.2.5流程图

开始

导入语音数据S

算法

数据归一化

建立滤波器重建语言

合成并加窗

确定帧数

线性预测法预测p个系数

计算

打开音频

男变女

开始

运行

女变男变童声

3数字变声器的软件实现及仿真结果

3.1 界面设计

在设计数字变声器的界面时，使用了MATLAB提供的可视化的界面环境Guide。相比通过编写程序进行GUI的设计，可视化的界面环境更加方便快捷。见图4

1按钮组：设计中使用了两个按钮组，把打开、原声、男变女、女变男、童声三项功能放在一组；另一个把语音的图像即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;

end

str=[pathname filename];

[temp,Fs]=audioread(str);

减小基音周期

提取共振峰

合成激励合成语音

图4 界面设计

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); %计算基音周期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 %增加共振峰

if imag(poles(p))>0

poles(p) = poles(p)*exp(1j*deltaOMG);

elseif imag(poles(p))<0

poles(p) = poles(p)*exp(-1j*deltaOMG);

end

end

end

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.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:143

R(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*Rop

Rop=R2;

Top=T2;

end

if R3>0.85*Rop

Rop=R3;

Top=T3;

end

PT=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); %计算基音周期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 %增加共振峰

if imag(poles(p))>0

poles(p) = poles(p)*exp(1j*deltaOMG);

elseif imag(poles(p))<0

poles(p) = poles(p)*exp(-1j*deltaOMG);

end

end

end

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

%绘图

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); %计算基音周期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/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

end

end

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.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）男声变童声图形对比

图5男声变女声前后图形对比

图6女声变男声前后图形对比

图6男声变童声前后图形对比

4 结束语

在开始制作时毫无头绪，但经过老师和同学的讲解，及自己所找到的各种文献及例子一点一点的探索，最终一点点做出来了。但还是有一点瑕疵，正努力找资料修改。在制作中遇到了很多比较棘手的问题如在滤波器参数的选取上，代码的更改问题，经常运行报错等。但最终经过同学和我自己的共同努力终于成功了。另外，在设计滤波器的过程中我们采用的是窗函数法，事实上，采用双线性变换法也可以进行，因此，不管对于多么简单的课题，其实也是有很多东西可以发掘的，这需要我们在平时多积累，多思考，只有这样，才能取得更大的进步，才能学有所用，学有所长。

参考文献

[1] 高西全、丁玉美编著.数字信号处理.西安.西安电子科技大学出版社，2008.

[2] 刘树棠译.数字信号处理——使用MATLAB.西安：西安交通大学出版社，2002.

[3] 胡广书编著.数字信号处理——理论、算法与实现.北京：清华大学出版社，2002.

[4] 梁虹等编.信号与线性系统分析————基于MATLAB的方法与实现.北京：高等教育出版社，2006.

[5] 郑阿奇.MATLAB实用教程。北京：电子工业出版社，2007.8

[6]蔡莲红，黄德智，蔡锐.现代语音技术基础与应用.北京：清华大学出版社,2003.

[7]张雄伟，陈亮，杨吉斌.现代语音处理技术及应用.北京：机械工业出版社，2003.