《MATLAB语音信号分析与合成（第二版）》：第3章 语音信号在其他变换域中的分析技术和特性

• 前言
• 1. 数据与函数路径设置
• 2. MATLAB仿真一：信号倒谱图、声道冲激响频谱、声门激励脉冲频谱
• 3. MATLAB仿真二：离散余弦逆变换
• 4. MATLAB仿真三：Mel滤波器组的频率响应曲线
• 5. MATLAB仿真四：MFCC参数匹配
• 6. MATLAB仿真五：小波分解
• 7. MATLAB仿真六：小波包分解
• 8. MATLAB仿真七：EMD分解
• 小结

前言

《MATLAB语音信号分析与合成（第二版）》是中科院声学所的大佬宋知用老师数十年经验积累下的呕心之作，对于语音信号处理相关感兴趣的同学，日后希望在语音信号分析、处理与合成相关领域进行一定研究的话，可以以此进行入门。

语音信号处理是数字信号处理的一个重要分支。本书含有许多数字信号处理的方法和 MATLAB函数。 全书共10章。第1-4章介绍语音信号处理的一些基本分析方法和手段,以及相应的MATLAB函数;第5-9章介绍语音信号预处理和特征的提取,包括消除趋势项和基本的减噪方法,以及端点检测、基音的提取和共 振峰的提取,并利用语音信号处理的基本方法,给出了多种提取方法和相应的 MATLAB程序;第10章结合 各种参数的检测介绍了语音信号的合成、语音信号的变速和变调处理,还介绍了时域基音同步叠加( TD PSOLA)的语音合成,并给出了相应的MATLAB程序。附录A中给出了调试复杂程序的方法和思路。 本书可作为从事语音信号处理的本科高年级学生、研究生或科研工程技术人员的辅助读物,也可作为从 事信号处理研究与应用的科研工程技术人员的参考用书。

1. 数据与函数路径设置

书中经常会调用的一些函数(自编函数或取自其他应用工具箱中的函数)已集中在basic_tbx工具箱中，在运行本书的程序前请把该工具箱设置(用set path设置)在工作路径下；

当要运行EMD处理时，要把emd工具箱设置在工作路径下;

当要运行主体延伸基音检测时，要把Pitch_ztlib工具箱设置在工作路径下;

当要进行时域基音同步叠加语音合成时，要把psola_lib工具箱设置在工作路径下;

当要应用本书提供的语音数据时，最好把speech_signal设置在工作路径下。

本书的所有函数和程序都在MATLAB R2009a版本下调试通过。（我用的是MATLAB2015b，有些函数已经更新，所以我会进行修改，以便调试通过）

路径设置的方法如下：

打开MATLAB，点击“主页”，找到设置路径

将上述文件夹路径全部添加到MATLAB搜索路径中

添加完毕，保存，开始仿真。

2. MATLAB仿真一：信号倒谱图、声道冲激响频谱、声门激励脉冲频谱

%% pr3_1_1clear all; clc; close all;y=load('su1.txt');% 读入数据fs=16000; nfft=1024;     % 采样频率和FFT的长度time=(0:nfft-1)/fs;      % 时间刻度figure(1), subplot 211; plot(time,y,'k');     % 画出信号波形title('信号波形'); axis([0 max(time) -0.7 0.7]);ylabel('幅值'); xlabel(['时间/s' 10 '(a)']); grid;figure(2)nn=1:nfft/2; ff=(nn-1)*fs/nfft; % 计算频率刻度Y=log(abs(fft(y)));      % 按式(3-1-8)取实数部分subplot 211; plot(ff,Y(nn),'k'); hold on;     % 画出信号的频谱图z=ifft(Y); % 按式(3-1-8)求取倒谱figure(1), subplot 212; plot(time,z,'k');     % 画出倒谱图title('信号倒谱图'); axis([0 time(512) -0.2 0.2]); grid; ylabel('幅值'); xlabel(['倒频率/s' 10 '(b)']);mcep=29;   % 分离声门激励脉冲和声道冲激响应zy=z(1:mcep+1);zy=[zy' zeros(1,nfft-2*mcep-1) zy(end:-1:2)']; % 构建声道冲激响应的倒谱序列ZY=fft(zy);% 计算声道冲激响应的频谱figure(2), % 画出声道冲激响应的频谱，用灰线表示line(ff,real(ZY(nn)),'color',[.6 .6 .6],'linewidth',3);grid; hold off; ylim([-4 5]);title('信号频谱（黑线）和声道冲激响频谱（灰线）')ylabel('幅值'); xlabel(['频率/Hz' 10 '(a)']); ft=[zeros(1,mcep+1) z(mcep+2:end-mcep)' zeros(1,mcep)]; % 构建声门激励脉冲的倒谱序列FT=fft(ft);      % 计算声门激励脉冲的频谱subplot 212; plot(ff,real(FT(nn)),'k'); grid;% 画出声门激励脉冲的频谱title('声门激励脉冲频谱')ylabel('幅值'); xlabel(['频率/Hz' 10 '(b)']); 

3. MATLAB仿真二：离散余弦逆变换

%% pr3_2_1clear all; clc; close all;f=50;   % 信号频率fs=1000;% 采样频率N=1000; % 样点总数n=0:N-1;xn=cos(2*pi*f*n/fs);  % 构成余弦序列y=dct(xn) ;    % 离散余弦变换num=find(abs(y)<5);   % 寻找余弦变换后幅值小于5的区间y(num)=0;      % 对幅值小于5的区间的幅值都置为0zn=idct(y);    % 离散余弦逆变换subplot 211; plot(n,xn,'k'); % 绘制xn的图title('x(n)'); xlabel(['样点' 10 '(a)']); ylabel('幅值');subplot 212; plot(n,zn,'k'); % 绘制zn的图title('z(n)'); xlabel(['样点' 10 '(b)']); ylabel('幅值');% 计算重建率rp=100-norm(xn-zn)/norm(xn)*100

rp =   84.3566

4. MATLAB仿真三：Mel滤波器组的频率响应曲线

% pr3_3_1 clear all; clc; close all;% 调用melbankm函数,在0-0.5区间设计24个Mel滤波器,用三角形窗函数bank=melbankm(24,256,8000,0,0.5,'t');bank=full(bank);bank=bank/max(bank(:));% 幅值归一化df=8000/256;    % 计算分辨率ff=(0:128)*df;  % 频率坐标刻度for k=1 : 24    % 绘制24个Mel滤波器响应曲线    plot(ff,bank(k,:),'k','linewidth',2); hold on;endhold off; grid;xlabel('频率/Hz'); ylabel('相对幅值')title('Mel滤波器组的频率响应曲线')

5. MATLAB仿真四：MFCC参数匹配

%% pr3_3_2 clear all; clc; close all;% [x1,fs]=wavread('s1.wav');      % 读入信号s1-\i1\% x2=wavread('s2.wav');    % 读入信号s2-\i2\% x3=wavread('a1.wav');    % 读入信号a1-\a1\[x1,fs]=audioread('s1.wav');      % 读入信号s1-\i1\x2=audioread('s2.wav');    % 读入信号s2-\i2\x3=audioread('a1.wav');    % 读入信号a1-\a1\wlen=200;  % 帧长inc=80;    % 帧移x1=x1/max(abs(x1));      % 幅值归一化x2=x2/max(abs(x2));x3=x3/max(abs(x3));% 计算/i1/与/i2/之间的匹配比较[Dcep,Ccep1,Ccep2]=mel_dist(x1,x2,fs,16,wlen,inc);figure(1)plot(Ccep1(3,:),Ccep2(3,:),'k+'); hold onplot(Ccep1(7,:),Ccep2(7,:),'kx'); plot(Ccep1(12,:),Ccep2(12,:),'k^');plot(Ccep1(16,:),Ccep2(16,:),'kh'); legend('第3帧','第7帧','第12帧','第16帧',2)% legend('第3帧','第7帧','第12帧','第16帧')xlabel('信号x1');ylabel('信号x2')axis([-12 12 -12 12]);line([-12 12],[-12 12],'color','k','linestyle','--');title('/i1/与/i2/之间的MFCC参数匹配比较')% 计算/i1/与/a1/之间的匹配比较[Dcep,Ccep1,Ccep2]=mel_dist(x1,x3,fs,16,wlen,inc);figure(2)plot(Ccep1(3,:),Ccep2(3,:),'k+'); hold onplot(Ccep1(7,:),Ccep2(7,:),'kx'); plot(Ccep1(12,:),Ccep2(12,:),'k^');plot(Ccep1(16,:),Ccep2(16,:),'kh'); legend('第3帧','第7帧','第12帧','第16帧',2)% legend('第3帧','第7帧','第12帧','第16帧')xlabel('信号x1');ylabel('信号x3')axis([-12 12 -12 12]);line([-12 12],[-12 12],'color','k','linestyle','--');title('/i1/与/a1/之间的MFCC参数匹配比较')

6. MATLAB仿真五：小波分解

%% pr3_4_1 clear all; clc; close all;% [x,fs]=wavread('awav.wav');      % 读入语音数据[x,fs]=audioread('awav.wav');      % 读入语音数据N=length(x);% 信号长度x=x-mean(x);% 消除直流分量J=2; % 设小波变换级数为J[C,L] = wavedec(x,J,'db1');      % 对时间序列进行一维多分辨分解CaLen=N/2.^J;      % 估计近似部分的系数长度Ca=C(1:CaLen);     % 取近似部分的系数Ca=(Ca-min(Ca))./(max(Ca)-min(Ca));     % 对近似部分系数做规正处理for i=1:CaLen      % 对近似部分系数做削波    if(Ca(i)<0.8), Ca(i)=0; endend[K,V]=findpeaks(Ca,[],6); % 寻找峰值位置和数值lk=length(K);if lk~=0    for i=2 : lk dis(i-1)=K(i)-K(i-1)+1;  % 寻找峰值之间的间隔    end    distance=mean(dis);   % 取间隔的平均值    pit=fs/2.^J/distance  % 计算这一帧的基音频率else    pit=0;end% 作图subplot 211; plot(x,'k'); title('一帧语音信号')subplot 212; plot(Ca,'k');title('用小波分解得到的近似系数中心削波后的峰值图')

pit =   100

7. MATLAB仿真六：小波包分解

%% pr3_4_2 clear all; clc; close all;filedir=[]; % 设置语音文件路径filename='aa.wav'; % 设置文件名fle=[filedir filename]    % 构成路径和文件名的字符串% [xx, fs, nbits]=wavread(fle);    % 读入语音文件[xx, fs]=audioread(fle);    % 读入语音文件x=xx-mean(xx);     % 消除直流分量x=x/max(abs(x));   % 幅值归一化N=length(x);% 取信号长度T=wpdec(x,5,'db2');% 对时间序列进行一维小波包分解% 按指定的结点,对时间序列分解的一维小波包系数重构y(1,:)=wprcoef(T,[5 0]);y(2,:)=wprcoef(T,[5 1]);y(3,:)=wprcoef(T,[5 2]);y(4,:)=wprcoef(T,[5 3]);y(5,:)=wprcoef(T,[5 4]);y(6,:)=wprcoef(T,[5 5]);y(7,:)=wprcoef(T,[5 6]);y(8,:)=wprcoef(T,[5 7]);y(9,:)=wprcoef(T,[4 4]);y(10,:)=wprcoef(T,[4 5]);y(11,:)=wprcoef(T,[5 11]);y(12,:)=wprcoef(T,[5 12]);y(13,:)=wprcoef(T,[4 7]);y(14,:)=wprcoef(T,[3 4]);y(15,:)=wprcoef(T,[3 5]);y(16,:)=wprcoef(T,[3 6]);y(17,:)=wprcoef(T,[3 7]);% 作图subplot 511; plot(x,'k');ylabel('/a/'); axis tightfor k=1 : 4    subplot(5,2,k*2+1); plot(y((k-1)*2+1,:),'k');    ylabel(['y' num2str((k-1)*2+1)]); axis tight;    subplot(5,2,(k+1)*2); plot(y(k*2,:),'k');    ylabel(['y' num2str(k*2)]); axis tight;endfigure    for k=1 : 4    subplot(5,2,(k-1)*2+1); plot(y((k-1)*2+9,:),'k');    ylabel(['y' num2str((k-1)*2+9)]); axis tight;    subplot(5,2,k*2); plot(y(k*2+8,:),'k');    ylabel(['y' num2str(k*2+8)]); axis tight;endsubplot(5,2,9); plot(y(17,:),'k');ylabel('y17'); axis tight

8. MATLAB仿真七：EMD分解

%% pr3_5_1 clear all; clc; close all;fs=5000;  % 采样频率N=500;    % 样点数n=1:N;t1=(n-1)/fs;     % 设置时间x1=sin(2*pi*50*t1);     % 产生笫1个正弦信号x2=(1/3)*sin(2*pi*150*t1);     % 产生笫2个正弦信号z=x1+x2;  % 把两个信号叠加imp=emd(z);      % 对叠加信号进行EMD分解[m,n]=size(imp); % 求取EMD分解成几个分量% 作图subplot(m+1,1,1);% 画叠加信号plot(t1,z,'k');title('原始信号'); ylabel('幅值')subplot 312;     % 画第1个正弦信号line(t1,x2,'color',[.6 .6 .6],'linewidth',5); hold onsubplot 313;     % 画第2个正弦信号line(t1,x1,'color',[.6 .6 .6],'linewidth',5); hold onfor i=1:m    subplot(m+1,1,i+1); % 画EMD分解后的信号    plot(t1,imp(i,:),'k','linewidth',1.5); ylabel('幅值')    title(['imf' num2str(i)]);endxlabel('时间/s');

小结

语音信号的处理主要就在变换域中的信号处理，通过将语音信号变换到其他域中进行分析，本章节中主要对信号倒谱图、声道冲激响频谱、声门激励脉冲频谱、离散余弦逆变换、Mel滤波器组的频率响应曲线、MFCC参数匹配、MFCC参数匹配、小波分解、小波包分解、EMD分解进行介绍。

MFCC参数是语音识别中一个非常重要的参数，自己研究生期间就是通过对比LPCC参数与MFCC参数进行语音信号的特征提取，在MFCC参数上的匹配效果确实优于LPCC参数，而且开源语音识别工具kaldi中也是对MFCC进行了优化，kaldi之父目前就在小米继续升级开发kaldi，其中有很多都已经应用到小米的语音识别场景中。对本章内容感兴趣或者想充分学习了解的，建议去研习书中第三章节的内容。后期会对其中一些知识点在自己理解的基础上进行讨论补充，欢迎大家一起学习交流。

关于宋老师：宋知用——默默传授MATLAB与信号处理知识的老人家

本系列文章列表如下：
《MATLAB语音信号分析与合成（第二版）》：第2章 语音信号的时域、频域特性和短时分析技术
《MATLAB语音信号分析与合成（第二版）》：第3章 语音信号在其他变换域中的分析技术和特性
《MATLAB语音信号分析与合成（第二版）》：第4章 语音信号的线性预测分析
《MATLAB语音信号分析与合成（第二版）》：第5章 带噪语音和预处理
《MATLAB语音信号分析与合成（第二版）》：第6章 语音端点的检测（1）
《MATLAB语音信号分析与合成（第二版）》：第6章 语音端点的检测（2）
《MATLAB语音信号分析与合成（第二版）》：第7章 语音信号的减噪
《MATLAB语音信号分析与合成（第二版）》：第8章 基音周期的估算方法
《MATLAB语音信号分析与合成（第二版）》：第9章 共振峰的估算方法
《MATLAB语音信号分析与合成（第二版）》：第10章 语音信号的合成算法