语音合成.rar_signals and systems_语音处理

function speechproc() % 定义常数 FL = 80; % 帧长 WL = 240; % 窗长 P = 10; % 预测系数个数 s = readspeech('voice.pcm',100000); % 载入语音s L = length(s);% 读入语音长度 wi= zeros(P,1); 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); zf_syn=zeros(P,1); % 变调不变速滤波器 exc_syn_t = zeros(L,1); % 合成的激励信号（脉冲串） s_syn_t = zeros(L,1); zi_syn_t=zeros(P,1); zf_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); zf_syn_v=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); %用线性预测法计算P个预测系数 % A是预测系数，E会被用来计算合成激励的能量 if n == 27 zplane(1,A); end s_f = s((n-1)*FL+1:n*FL); % 本帧语音，下面就要对它做处理 % (4) 在此位置写程序，用filter函数s_f计算激励，注意保持滤波器状态 [w,wf]=filter(A,1,s_f,wi); exc((n-1)*FL+1:n*FL) =w; %滤波器的下一初态为当前末态 wi=wf; % (5) 在此位置写程序，用filter函数和exc重建语音，注意保持滤波器状态 [ww, zi_rec]=filter(1,A,exc((n-1)*FL+1:n*FL),zi_pre); s_rec((n-1)*FL+1:n*FL)=ww; zi_pre=zi_rec; % s_rec((n-1)*FL+1:n*FL) = ... 将你计算得到的重建语音写在这里 % 注意下面只有在得到exc后才会计算正确 s_Pitch = exc(n*FL-222:n*FL); PT = findpitch(s_Pitch); % 计算基音周期PT（不要求掌握） G = sqrt(E*PT); % 计算合成激励的能量G（不要求掌握） % (10) 在此位置写程序，生成合成激励，并用激励和filter函数产生合成语音 for i=(n-1)*FL+1:n*FL exc_syn(i)=G*(mod(i,PT)==0); end [s_syn1, zf_syn]=filter(1,A,exc_syn((n-1)*FL+1:FL*n),zi_syn); zi_syn=zf_syn; s_syn((n-1)*FL+1:n*FL)=s_syn1; % exc_syn((n-1)*FL+1:n*FL) = ... 将你计算得到的合成激励写在这里 % s_syn((n-1)*FL+1:n*FL) = ... 将你计算得到的合成语音写在这里 % (11) 不改变基音周期和预测系数，将合成激励的长度增加一倍，再作为filter % 的输入得到新的合成语音，听一听是不是速度变慢了，但音调没有变。 for i=(n-1)*FL*2+1:n*FL*2 exc_syn_v(i)=G*(mod(i,PT)==0); end [s_syn1, zf_syn_v]=filter(1,A,exc_syn_v((n-1)*FL*2+1:FL*n*2),zi_syn_v); zi_syn_v=zf_syn_v; s_syn_v((n-1)*FL*2+1:n*FL*2)=s_syn1; % exc_syn_v((n-1)*FL_v+1:n*FL_v) = ... 将你计算得到的加长合成激励写在这里 % s_syn_v((n-1)*FL_v+1:n*FL_v) = ... 将你计算得到的加长合成语音写在这里 % (13) 将基音周期减小一半，将共振峰频率增加150Hz，重新合成语音，听听是啥感受～ temp=floor(PT/2); for i=(n-1)*FL+1:n*FL exc_syn_t(i)=G*(mod(i,temp)==0); end [z,p,k]=tf2zp(1,A);%通过该函数计算极点位置 for i=1:length(p) if(angle(p(i))>0) p(i)=p(i)*exp(1i*150*2*pi/8000); else p(i)=p(i)*exp(-1i*150*2*pi/8000); end end [b,a]=zp2tf(z,p,k); [s_syn1,zf_syn_t]=filter(b,a,exc_syn_t((n-1)*FL+1:n*FL),zi_syn_t); zi_syn_t=zf_syn_t; s_syn_t((n-1)*FL+1:n*FL)=s_syn1; % exc_syn_t((n-1)*FL+1:n*FL) = ... 将你计算得到的变调合成激励写在这里 % s_syn_t((n-1)*FL+1:n*FL) = ... 将你计算得到的变调合成语音写在这里 end % (6) 在此位置写程序，听一听 s ，exc 和 s_rec 有何区别，解释这种区别 % 后面听语音的题目也都可以在这里写，不再做特别注明 sound(s/max(s),8000); pause(2); sound(exc/max(exc),8000); pause(2); sound(s_rec/max(s_rec),8000); pause(2) sound(s_syn/max(s_syn),8000); pause(2) sound(s_syn_v/max(s_syn_v),8000); pause(4) sound(s_syn_t/max(s_syn_t),8000); figure subplot(3,2,1),plot(s),title('s'); subplot(3,2,2),plot(exc),title('exc'); subplot(3,2,3),plot(s_rec),title('s_rec'); subplot(3,2,4),plot(s_syn),title('s_syn'); subplot(3,2,5),plot(s_syn_v),title('s_syn_v'); subplot(3,2,6),plot(s_syn_t),title('s_syn_t'); % 保存所有文件 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 % 从PCM文件中读入语音 function s = readspeech(filename, L) fid = fopen(filename, 'r'); s = fread(fid, L, 'int16'); fclose(fid); return % 写语音到PCM文件中 function writespeech(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