matlab高级编程与应用-语音处理实验报告

自03 张驰昱2010012028

一、 语音预测模型

(1) 给定e(n) = s(n) -a1s(n -1) -a2s(n -2)假设e(n)是输入信号，s(n)是输出信号，上述滤

波器的传递函数是什么？如果a1 = 1.37，a2 = -0.9506 ，上述合成模型的共振峰频率是多少？用zplane,freqz,impz分别绘出零极点图，频率响应和单位样值响应。用filter绘出单位样值响应，比较和impz的是否相同。

问题分析：

本问题主要练习传递函数到零极点的转化，零极点的绘制，频率响应的绘制，单位响应的绘制，复习filter数字滤波器的使用。

具体实现： clear;clc;

a = [1, -1.37, 0.9506];

sys=tf(1,a,-1,'variable','z^-1') [z,p]=tf2zp(1,a);

%[r,p,k]=residuez(1,a);也能求出零点 omg=abs(angle(p(1)));

fs=8000;%数字采样频率

f=omg*fs/2/pi%弧度转化为频率 n=[0:49]';x=(n==0); figure(1);zplane(1,a); figure(2);freqz(1,a);

figure(3);subplot(2,1,1),stem(n,filter(1,a,x)); figure(3);subplot(2,1,2),impz(1,a,50);

(2) 理解speechproc的主要流程

我认为主要的部分是以下程序段：(个人理解写在了注释中) %先要统一初始化所用到的向量，这样可以提高执行效率

for n = 3:FN%汉明窗取到了帧长的三倍，所以n从3开始 s_w = s(n*FL-WL+1:n*FL).*hw; %加窗方便用lpc处理

[A E] = lpc(s_w, P); %用lpc技术得到传递函数系数A s_f = s((n-1)*FL+1:n*FL); %待处理的本帧语音，即激励响应 %需要推算本帧语音的激励，只有得到了激励才能做接下来的变声处理

s_Pitch = exc(n*FL-222:n*FL); PT= findpitch(s_Pitch); %刚才算出的激励信号是有高斯白噪声的，需要找

%出基音周期和能量，为重新合成激励信号做准备

G = sqrt(E*PT(n));

(3) 在27帧处观察零极点图

问题分析：主要让我们对语音传函的共轭极点有一个更直观的认识

具体实现： if n == 27

figure(n);zplane(1,A); end

(4)用filter计算每帧的激励信号

问题分析：已经求出了传函系数和激励相应，只要传函的分子分母互换把激励相应当激

励，得到的相应就是原激励

具体实现：

%前输出状态作为后输入状态即前后状态不变

[temp1,zi_pre]=filter(A,1,s_f,zi_pre); exc((n-1)*FL+1:n*FL)=temp1;

(5)利用刚才得到的激励信号，继续用filter重建语音

问题分析：相当于对于之前求出的激励的验算。与前题一样，要保持前后状态不变。

具体实现：

[temp2,zi_rec]=filter(1,A,temp1,zi_rec); s_rec((n-1)*FL+1:n*FL)=temp2; (6)比较激励信号与s(n)和重建语音的区别

问题分析：重建语音因为用原激励exc做输入，所以和原语音相同。但exc中包含了高斯白噪声序列，这就是为什么接下来要计算exc的基音周期并重新合成“干净”的激励。

具体实现：

figure;

subplot(2,1,1),plot(s(2000:4000)); subplot(2,1,2),plot(exc(2000:4000));

二、 语音合成模型

(1)

生成一个8kHz抽样频率的持续1s的数字信号，该信号是频率为200Hz的单位样值串

问题分析：200Hz即一秒钟200个样值，所以NS=200。又因为抽样频率为8000，所以相邻样值间隔为40，即N=40。

具体实现：

N=8000;

f=200; %如果要300Hz，就改为300 e=(mod([1:N],floor(N/f))==0);

%8000个点取模再求逻辑判断自然生成一个样值串 sound(e,8000);

(2)

生成基因周期满足PT=80+5mod(m,50)的样值序列，其中m为10ms长的段的序号。

问题分析：需要有一个变量来记忆现在已经到了第几个段，用循环累加的方式存储之前所有样值间隔的和。而且没法用for循环，因为步距不固定且步数事先不知道。用while循环。

具体实现： N=8000;

e=zeros(N,1); %一秒信号初始化 PTsum=1; %第一个样值在1处 whilePTsum<=N

e(PTsum)=1;

m=ceil(PTsum/80); %当前位置所处的段序号，注意这里是向上取整 PT=80+5*mod(m,50); %当前段的基音周期 PTsum=PTsum+PT; %累加样值间隔 end

sound(e,8000);

(3)

把刚生成的激励放到第一题的系统中。试听激励与响应的区别

具体实现：

在上题程序后加上一下语句即可 a = [1, -1.37, 0.9506]; s=filter(1,a,e); sound(s,8000);

感觉激励是憋在喉咙里的声音，这里的响应是张开嘴的声音

(4)

用每一帧得到的基因周期和(2)的方法，合成激励Gx(n),再用filter还原语音。比较与原语音的差别。

问题分析：仍用PTsum来记录当前位置，由于PT已求得，所以比(2)简单一些。注意PTsum的初始化要放在大循环外。

具体实现：

PTsum=2*FL+1; %大循环从第3帧开始，所以第一个样值是2FL+1点 for n = 3:FN

„„

whilePTsum<=n*FL exc_syn(PTsum)=G;

PTsum=PTsum+PT;

End%样值串生成完毕

[temp3,zi_syn]=filter(1,A,exc_syn((n-1)*FL+1:n*FL),zi_syn); s_syn((n-1)*FL+1:n*FL)=temp3; %本帧重建语音 „„ end

与原语音相比，此次重建的语音噪声更小。

三、 变速不变调

把每帧激励增长一倍，生成比原语音慢一倍的语音。

问题分析：

只要修改生成激励的帧长即可。和上题基本一样。

具体实现：

FL_v=FL*2;%在循环外定义新的FL_v和PTsum_v PTsum_v=2*FL_v+1; for n = 3:FN

„„

whilePTsum_v<=n*FL_v

exc_syn_v(PTsum_v)=G; %注意激励的幅度为G PTsum_v=PTsum_v+PT; end

[temp3,zi_syn_v]=filter(1,A,exc_syn_v((n-1)*FL_v+1:n*FL_v),zi_syn_v);

s_syn_v((n-1)*FL_v+1:n*FL_v)=temp3; „„ end

四、变速不变调 (1)

将第一题中的系统的共振峰频率提高150Hz后a1和a2是多少。

问题分析：

主要练习传函到零极点，再从零极点到传函的转化。注意是极点幅角的绝对值增大，

横轴上的极点不动。

具体实现： f_plus=150;

a = [1 -1.37 0.9506];

[z,p,k]=tf2zp(1,a); %传函到零极点

p=p.*exp(j*f_plus*pi/4000*sign(angle(p)).*(angle(p)代码优化：

一开始认为要用循环加判断来旋转极点，但是思考后发觉这也是可以转化为矩阵的： p=p.*exp(j*f_plus*pi/4000*sign(angle(p)).*(angle(p)(2)

将基音周期减小一半，共振峰频率增加150Hz

问题分析：

在重建语音之前先重新计算传函系数A，这在上问已经实现；在合成激励时PT减小一般即可。注意用来记忆当前位置的PTsum要用新的，不可沿用上题的。即重建几遍语音就有几个不相关的PTsum

具体实现：

f_plus=150;%共振峰提高频率 PTsum_t=2*FL+1;

times=2;%基因周期减小倍数 for n = 3:FN

„„

[z,p,k]=tf2zp(1,A);

p=p.*exp(sign(angle(p)).*(angle(p)PTsum_t=PTsum_t+floor(PT/times); %注意要取整 End

[temp3,zi_syn_t]=filter(1,A,exc_syn_t((n-1)*FL+1:n*FL),zi_syn_t);

s_syn_t((n-1)*FL+1:n*FL)=temp3; „„ end

五、 逆向工程

（1） 已知有一款商业娱乐产品可实现语音信号的变调处理。现进行如下实验：将该产品

放置于音箱前，用计算机播放voice.pcm，声音经音箱放出后被该产品采集到，随后该产品生成变调语音、播放并自行记录，最后将其记录的音频文件导入计算机得到Tomvoice.pcm。请你在分析比对这两个语音信号的频谱特征的基础上，研究该产品生成变调语音的机理，画出实现框图，并指明具体参数。

解决步骤：

1、 选取两语音中说的同一个字分析其频谱，比较共振峰频率提高了多少。

2、 沿用speechproc中的程序段，分别计算两语音各帧基因周期的平均值。比较

处理前后基因周期相差多少。

3、 比较两信号长度得处理后提速多少。

具体实现： 步骤一：

l=200; %取25ms的语音样本

fid = fopen('voice.pcm', 'r'); s = fread(fid,inf,'int16'); fclose(fid);

f(:,1)=s(2428:2428+l-1); %选取了“电灯”的“灯”的前25ms fid = fopen('Tomvoice.pcm', 'r'); s_tom = fread(fid,inf,'int16'); fclose(fid)

f(:,2)=s_tom(2335:2335+l-1); %同样选取了“电灯”的“灯”的前25ms OMG=5000*pi;%频域宽度 k=1000;%频域采样点个数

omg=linspace(0,OMG-OMG/k,k);

t=linspace(0,(l-1)/8000,l);%时域宽度即音频播放时间

F=(l-1)/8000/l*exp(-j*kron(omg',t))*f; %两个信号一起做变换 F=abs(F); figure;

subplot(2,1,1),plot(omg/2/pi,F(:,1)); %横坐标为频率 subplot(2,1,2),plot(omg/2/pi,F(:,2));

运行结果：

从图上得峰值左移了263Hz

步骤二：

%沿用speechproc中的程序段

fid = fopen('filename.pcm', 'r'); s = fread(fid,inf,'int16'); fclose(fid); FL = 80; WL = 240; P = 10;

L = length(s);

FN = floor(L/FL)-2;

exc = zeros(L,1);

zi_pre = zeros(P,1); hw = hamming(WL); PT=zeros(FN,1); for n = 3:FN

s_w = s(n*FL-WL+1:n*FL).*hw;

[A E] = lpc(s_w, P);

[exc((n-1)*FL+1:n*FL),zi_pre]=filter(A,1,s((n-1)*FL+1:n*FL),zi_pre);

PT(n) = findpitch(exc(n*FL-222:n*FL)); %把每个PT都存在列向量里

End

PT_ava=sum(PT)/(FN-2)%计算PT平均值

分别对voice.pcm和Tomvoice.pcm运行后得到PT_ava分别为50.79和38.411，即基因周期缩短了1.32倍。 步骤三：

通过比较信号长度知，处理后语音提速了约1.1倍。

整体实现框图：

录音 提取录音信号 (3)

变声处理： ·提速1.1倍 ·共振峰频率提高263Hz ·基因周期缩短1.32倍 保存新语音 自制“趣味学舌”软件，可自动录音，也可手动录音。

问题分析：继续沿用speechproc函数。此问题即在变调又变速的基础上再合成录音的功能，所以主要的任务其实在编写录音函数。用一个参量来控制是自动还是手动录音。自动录音采用短时的循环来甄别是否振幅超过某标杆值，一旦侦测到，即开始录音。

具体实现：

functiontom()

auto=1;

%当auto为1时自动录音，即运行程序后，一有语音就开始录音。

%当auto为0时手动录音，即运行程序1s后开始录音。具体到录音函数中还有说明。 clc;clear;

FL = 80; WL = 240; P = 10; s = TomR(auto); L = length(s); FN = floor(L/FL)-2;

exc = zeros(L,1); %预测激励 exc_tom = zeros(L,1);%合成激励 s_tom = zeros(L,1); %新生成的tom语音 zi_pre = zeros(P,1);%原语音状态 zi_tom = zeros(P,1);%tom语音的状态 hw = hamming(WL);

speedup=1.1; %速度加快倍数 f_plus=263;%共振峰频率提高

short=1.32; %基因周期缩短倍数 FL_v=floor(FL/speedup);%新合成的帧长 PTsum=2*FL_v+1;%记录当前合成激励位置 for n = 3:FN

s_w = s(n*FL-WL+1:n*FL).*hw; [A E] = lpc(s_w, P);

[exc((n-1)*FL+1:n*FL),zi_pre]=filter(A,1,s((n-1)*FL+1:n*FL),zi_pre); PT= findpitch(exc(n*FL-222:n*FL)); G = sqrt(E*PT(n));

%语音变换环节现在开始，其实就是把变调变速合并起来 [z,p,k]=tf2zp(1,A);

p=p.*exp(sign(angle(p)).*(angle(p)PTsum=PTsum+floor(PT/short); end

[temp,zi_tom]=filter(1,A,exc_tom((n-1)*FL_v+1:n*FL_v),zi_tom); s_tom((n-1)*FL_v+1:n*FL_v)=temp; end

displayTom: sound(s_tom,8000);

fid = fopen('tom.pcm','w'); fwrite(fid, s_tom, 'int16'); fclose(fid); end

%录音函数

function s=TomR(auto)

t=3;%不管手动模式还是自动模式，开始录音后3秒钟停 if(auto==1)%自动模式 a=0;

r = audiorecorder(8000, 16, 1); record(r);

pause(0.5);%停0.5秒是因为开始刚录音时有个冲击的噪声 while(a<0.003)%循环检测有无语音，一旦检测到就跳出循环 r = audiorecorder(8000, 16, 1); record(r);

pause(0.01);%循环周期是0.01秒

pause(r);%只有暂停后才可获取语音，所以获得录音总会有一个空白间隙 s = getaudiodata(r, 'double');

a=mean(abs(s));%计算这0.01秒内振幅的平均值 end

displayStart!

resume(r);%为了尽可能的保住语音，用resume而不重新生成 pause(t);%录三秒 stop(r); displayDone!

s = getaudiodata(r, 'double'); end

if(auto==0)%手动模式

r = audiorecorder(8000, 16, 1); displayReady

pause(1);%准备一秒再开始 displayStart! record(r);

pause(t);%录三秒 displayDone. stop(r);

s = getaudiodata(r, 'double');%获取录音信号 end end

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; end

设auto=1后，除了看不见那只猫以外，感觉跟iphone上的Talking Tom很像。当调试完这个程序，（主要在调TomR这个录音函数），天已经亮了。但是当我对着屏幕说话，matlab会自动开始录音并回话的时候，感觉心情大爽。

六、 总结

此次试验所需要的写的程序不是太长，但是我觉得新学到的知识很丰富。首先让我对发声机理和语音处理有一个比较入门和宏观的认识，虽然很关键的lpc技术和基音频率寻找算法没有掌握，但是用一个离散序列的传递函数和一个样值序列来描述声门脉冲通过声管的过程让我感到学以致用的快感。还了解到不同人的声音的区别其实就是不同声管传递函数极点的位置区别和声门脉冲基音周期的区别。对离散信号的传递理解的更直观了。

附：完整的speechproc函数，相关部分的注释已在前文全部解释过。

functionspeechproc() clc;clear;

FL = 80; WL = 240; P = 10;

s = readspeech('voice.pcm',inf); L = length(s); FN = floor(L/FL)-2;

exc = zeros(L,1); zi_pre = zeros(P,1);

s_rec = zeros(L,1); zi_rec = zeros(P,1); exc_syn = zeros(L,1); s_syn = zeros(L,1); zi_syn = zeros(P,1); exc_syn_t = zeros(L,1); s_syn_t = zeros(L,1); zi_syn_t=zeros(P,1); exc_syn_v = zeros(2*L,1); s_syn_v = zeros(2*L,1); zi_syn_v=zeros(P,1); hw = hamming(WL); PTsum=2*FL+1; FL_v=FL*2; PTsum_v=2*FL_v+1;

PT=zeros(FN,1); f_plus=150; PTsum_t=2*FL+1;

times=2; for n = 3:FN

s_w = s(n*FL-WL+1:n*FL).*hw; [A E] = lpc(s_w, P); if n == 27 %(3)

figure(n);zplane(1,A); end

s_f = s((n-1)*FL+1:n*FL); %(4)

[temp1,zi_pre]=filter(A,1,s_f,zi_pre); exc((n-1)*FL+1:n*FL)=temp1; %(5)

[temp2,zi_rec]=filter(1,A,temp1,zi_rec); s_rec((n-1)*FL+1:n*FL)=temp2;

s_Pitch = exc(n*FL-222:n*FL); PT(n) = findpitch(s_Pitch); G = sqrt(E*PT(n)); %(10)

whilePTsum<=n*FL exc_syn(PTsum)=G; PTsum=PTsum+PT(n); end

[temp3,zi_syn]=filter(1,A,exc_syn((n-1)*FL+1:n*FL),zi_syn); s_syn((n-1)*FL+1:n*FL)=temp3; % (11)

whilePTsum_v<=n*FL_v exc_syn_v(PTsum_v)=G; PTsum_v=PTsum_v+PT(n); end

[temp3,zi_syn_v]=filter(1,A,exc_syn_v((n-1)*FL_v+1:n*FL_v),zi_syn_v); s_syn_v((n-1)*FL_v+1:n*FL_v)=temp3;

%(13)

[z,p,k]=tf2zp(1,A);

p=p.*exp(sign(angle(p)).*(angle(p)PTsum_t=PTsum_t+floor(PT(n)/times); end

[temp3,zi_syn_t]=filter(1,A,exc_syn_t((n-1)*FL+1:n*FL),zi_syn_t); s_syn_t((n-1)*FL+1:n*FL)=temp3; end %(6)

sound(s,8000); sound(exc,8000); sound(s_rec,8000); sound(s_syn,8000); sound(s_syn_v,8000); sound(s_syn_t,8000); figure;

subplot(2,1,1),plot(s(2000:4000)); subplot(2,1,2),plot(exc(2000:4000)); % ±£´æËùÓÐÎļþ

writespeech('exc.pcm',exc); writespeech('rec.pcm',s_rec); writespeech('exc_syn.pcm',exc_syn); writespeech('syn.pcm',s_syn);

writespeech('exc_syn_t.pcm',exc_syn_t); writespeech('syn_t.pcm',s_syn_t); writespeech('exc_syn_v.pcm',exc_syn_v); writespeech('syn_v.pcm',s_syn_v); return

function s = readspeech(filename, L) fid = fopen(filename, 'r'); s = fread(fid, L, 'int16'); fclose(fid); return

functionwritespeech(filename,s) fid = fopen(filename,'w'); fwrite(fid, s, 'int16'); fclose(fid); return

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