数字音效处理器 项目报告

1、数字信号处理B课程项目终期报告题 目： 数字音效处理器 组 号： 14 号 任课教师： 组 长： 成 员： 成 员： 成 员： 成 员： 联系方式： 二0一五 年 10月 24日目录一 项目介绍二 项目原理三 项目完成过程四 项目结果与分析五 工作分配六 参考文献一 项目介绍本次项目要求设计一个数字音效处理器，通过一系列相关算法以及借助相关工具（MATLAB），对声音信号进行处理，要求具备语音识别以及降除噪声的功能，另外我们又加入了语音均衡器，丰富了数字音效处理器的功能。最后加入GUI界面设计，方便了用户的使用。二 项目原理本部分主要分为三部分，分别是语音识别，噪声降除以及语音均衡器的相关原理

2、。语音识别：1.1项目大体步骤我们主要采用基于VQ的说话人识别，我们采用的识别结构框图如下，提取特征函数预处理语音信号VQ码本匹配计算VQ码本训练码本库识别结果1.2 语音信号的预处理预处理包括：预加重，端点检测，分帧以及加窗。预加重的目的是将更为有用的高频部分频谱进行提升，使信号的频谱变得平坦，保持在低频到高频的整个频带中，能用同样的信噪比求取频谱，以便于进行频谱分析或声道参数分析。端点检测即是对输入语音信号的起始点与结束点的判定。分帧：由于语音信号的准平稳特性，使得其只有在短时段上才可被视为一个平稳过程，所以需要把一定长度的语音分为很多帧来分析。加窗：为了减少语音帧的截断效应，降低帧两端的

3、坡度，使语音帧的两端不引起急剧变化而平滑过渡到零，需要将语音帧乘以一个窗函数。1.3 MFCC特征函数的提取Mel频率倒谱系数（MFCC）的分析与传统的线性倒谱系数（LPCC）不同，它的分析是着眼于人耳的听觉机理，因为人类在对1000Hz以下的声音频率范围的感知遵循近似的线性关系；对1000Hz以上的声音频率范围的感知不遵循线性关系，而是遵循在对数频率坐标上的近似线性关系，所以Mel 倒谱系数获得了较高的识别率和较好的鲁棒性。实现上，Mel倒谱系数是将语音频率划分成一系列三角形的滤波器序列，这组滤波器在频率的Mel坐标上是等待宽的。MFCC参数的提取过程为：1. 对输入的语音信号进行

4、分帧、加窗，然后作离散傅立叶变换，获得频谱分布信息。 设语音信号的DFT为： *（1）2. 再求频谱幅度的平方，得到能量谱。 3. 将能量谱通过一组Mel尺度的三角形滤波器组。 我们定义一个有M个滤波器的滤波器组（滤波器的个数和临界带的个数相近），采用的滤波器为三角滤波器，中心频率为f(m),m=1,2,3,···,M，本系统取M=100。4. 计算每个滤波器组输出的对数能量。 *（2）其中为三角滤波器的频率响应。 5. 经过离散弦变换（DCT）得到MFCC系数。*（3） MFCC系数个数通常取2030，常常不用0阶倒谱系数，因为它反映的是频谱能量，故在一般识别系

5、统中，将称为能量系数，并不作为倒谱系数，本文选取20阶倒谱系数。与线性预测倒谱分析相比，MFCC的优点是不依赖全极点语音信号产生模型的假定，因而在噪声环境下表现出更强的鲁棒性，在非特定人说话人识别方面有利于减小因说话人不同的差异可能带来的影响。1.4矢量量化矢量量化是将若干个幅度连续取值的时域采样信号分成一组，则构成矢量，然后用若干离散数字值，来表示各种矢量，其有效地利用了矢量中各分量的四种相关的性质，因此在量化过程中最大限度地消除数据之间的冗余，从而实现高效的熵压缩码。其基本原理为：将一帧语音数据中提取的矢量特征在多维空间给予整体量化。首先把信号序列的每隔K个连续点分成一组，形成K维欧氏空间

6、中的一个矢量，然后对此矢量进行量化。本文中矢量量化有两个作用，在训练阶段，把每个说话者所提取的特征参数进行分类，产生不同码子组成的码本；在识别阶段，我们用VQ方法计算平均失真测度，以此进行说话人识别。VQ码本的设计主要采用LBG算法，该算法是一种递推算法，从一个事先选定的初始码本进行迭代，直到系统性能满足要求或不再有明显的改进为止。具体LBG算法如下：设 为训练序列，B为码本。 1. 取提取出来的所有帧的特征矢量的型心(均值)作为第一个码字矢量B1。2. 将当前的码本Bm根据以下规则分裂，形成2m个码字。 *（4）其中m从1变化到当前的码本的码字数，是分裂时的参数，本文=0.01。 3. 根据

7、得到的码本把所有的训练序列(特征矢量)进行分类，然后按照下面两个公式计算训练矢量量化失真量的总和以及相对失真(n为迭代次数，初始n=0，=，B为当前的码书)，若相对失真小于某一阈值，迭代结束，当前的码书就是设计好的2m个码字的码书，转。否则，转下一步。 量化失真量和： 相对失真： *（6）4. 重新计算各个区域的新型心，得到新的码书，转3。 5. 重复 ， 和步，直到形成有M个码字的码书(M是所要求的码字数)，其中D0=10000。1.5说话人识别在匹配阶段，VQ用来计算未知说话人特征矢量和已经训练说话人码本之间的平均量化失真D，并设置一个阈值，若D小于此阈值，则是原训练者，反之则认为不是原训

8、练者。注： *（7）2噪声降除我们使用了维纳法和LMS自适应滤波法来降除噪声。2.1 维纳法（Wiener Scalart Method）维纳是用来解决从噪声中提取信号的一种过滤（或滤波）方法。这种线性滤波问题，可以看做是一种估计问题或一种线性估计问题。 一个线性系统，如果它的单位样本响应为h(n)，当输入一个随机信号x(n) 且x(n)=s(n)+v(n) *(8)其中x(n)表示信号，v(n)表示噪声，则输出y(n)为 yn=mh(m)x(n-m) *(9)我们希望x(n)通过线性系统h(n)后得到的y(n)尽量接近于s(n)，因此称y(n)为s(n)的估计值，用s(n)表示，即 yn=s

9、(n). *(10) 则维纳滤波器的输入输出关系可用下面图1表示。 实际上，式 (9) 所示的卷积形式可以理解为从当前和过去的观察值x(n)，x(n-1)，X(n-2)x(n-m).来估计信号的当前值 s(n)。因此，用h(n) 进行过滤问题实际上是一种统计估计问题。 一般地，从当前的和过去的观察值x(n)，x(n-1)，x(n-2)估计当前的信号值Y(n)=s(n)成为过滤或滤波；从过去的观察值，估计当前的或者将来的信号值 yn=s(n+N).称为外推或预测；从过去的观察值，估计过去的信号值 yn=s(n-N) 称为平滑或内插。因此维纳滤波器又常常被称为最佳线性过滤与预测或线性最优估计。这里

10、所谓的最佳与最优是以最小均方误差为准则的。 如果我们分别以s(n) 与 s(n) 表示信号的真实值与估计值，而用e(n) 表示他们之间的误差，即 en=sn-s(n) *(11)显然e(n)可能是正值，也可能是负值，并且它是一个随机变量。因此，用它的均方误差来表达误差是合理的，所谓均方误差最小即它的平方的统计期望最小： n=Ee2n=min* (12)采用最小均方误差准则作为最佳过滤准则的原因还在于它的理论分析比较简单，不要求对概率的描述。 为了按 (12) 式所示的最小均方误差准则来确定维纳滤波器的冲激响应h(n)，令n对h(j)的导数等于零，即可得 Rxsm=hiRxxm-1,m * (1

11、3)式中，Rxsm 是s(n)与x(n)的互相关函数，Rxxm是x(n)的自相关函数，分别定义为 Rxs=Exxsn+m Rxx=Ex(x)x(n+m)式 (13) 称为维纳滤波器的标准方程或维纳霍夫方程。如果已知 Rxsm 和 Rxxm,那么解此方程即可求的维纳滤波器的冲激响应。式 (13) 所示标准方程右端的求和范围即i的取值范围没有具体标明，实际上有三种情况：(1) 有限冲激响应（FIR）维纳滤波器，i从0到N-1取得有限个整数值；(2) 非因果无限冲激响应（非因果IIR）维纳滤波器，i从到取所有整数值；(3) 因果无限冲激响应（因果IIR）维纳滤波器，i从0到取正整数值。上述三种情况下

12、标准方程的解法不同，本文只描述FIR维纳滤波器的求解。设滤波器冲激响应序列的长度为N，冲激响应矢量为 h=h0h1hN-1T * (14)滤波器输入数据矢量为 xn=xnxn-1x(n-N+1)T *(15)则滤波器的输出为 yn=sn=xTnh=hTx(n) *(16)这样，式 (13) 所示的维纳-霍夫方程可写成 PT=hTR or P=Rh * (17)其中 P=Ex(n)s(n) * (18)是s(n) 与x(n) 的互相关函数，它是一个N维列矢量；R是x(n)的自相关函数，是N阶方阵R=Ex(n)xT(n) * (19)利用求逆矩阵的方法直接求解式 (19)，得 hopt=R-1P

13、*(20)这里opt表示“最佳”，这就是FIR维纳滤波器的冲激响应。2.2 LMS自适应滤波法自适应滤波器是能够根据输入信号自动调整性能进行数字信号处理的数字滤波器。作为对比，非自适应滤波器有静态的滤波器系数，这些静态系数一起组成传递函数。对于一些应用来说，由于事先并不知道所需要进行操作的参数，例如一些噪声信号的特性，所以要求使用自适应的系数进行处理。在这种情况下，通常使用自适应滤波器，自适应滤波器使用反馈来调整滤波器系数以及频率响应。总的来说，自适应的过程涉及到将代价函数用于确定如何更改滤波器系数从而减小下一次迭代过程成本的算法。价值函数是滤波器最佳性能的判断准则，比如减小输入信号中的噪声成

14、分的能力。随着数字信号处理器性能的增强，自适应滤波器的应用越来越常见，时至今日它们已经广泛地用于手机以及其它通信设备、数码录像机和数码照相机以及医疗监测设备中。下面图示的框图是最小均方滤波器（LMS）和递归最小平方（Recursive least squares filter，RLS，即我们平时说的最小二乘法) 这些特殊自适应滤波器实现的基础。框图的理论基础是可变滤波器能够得到所要信号的估计。1.输入信号是所要信号  d(n) 和干扰噪声  v(n) 之和: x(n)=d(n)+v(n) 2.可变滤波器有有限脉冲响应结构，这样结构的脉冲响应等于滤波器系数

15、。p 阶滤波器的系数定义为 Wn=n0,n1,npT *(21)3.误差信号或者叫作代价函数，是所要信号与估计信号之差en=dn-d(n) *(22)可变滤波器通过将输入信号与脉冲响应作卷积估计所要信号，用向量表示为 dn=wn*x(n) *(23)其中x(n)=x(n),x(n-1),x(n-p)T是输入信号向量。另外，可变滤波器每次都会马上改变滤波器系数 wn+1=wn+wn *(24)其中  wn是滤波器系数的校正因子。自适应算法根据输入信号与误差信号生成这个校正因子，LMS 和 RLS 是两种不同的系数更新算法。3语音均衡器我们所设计的数字均衡器包括以下功

16、能：1 音频文件（.wav）的获取；2 八段数字均衡器设定；3 滤波；4 播放，播放均衡过后的语音信号，并显示相关波形。此设计的关键是设计出满足一定性能指标的各个数字滤波器，构成一个覆盖音频频段的滤波器组。3.1数字滤波器数字滤波器的功能是把输入序列通过一定的运算，变换成输出序列。数字滤波器分为两类：无限冲激响应（IIR）以及有限冲激响应（FIR）（1） IIRIIR滤波器的最大特点是：单位冲激响应h（n）是无限长的；系统函数H（z）在有限长Z平面内有极点存在；结构上存在输出到输入的反馈，即结构上是递归的；稳定的IIR滤波器其全部极点一定在单位圆内。（2） FIRFIR滤波器最大的特点是：系统

17、的单位冲激响应h（n）是个有限长序列；系统函数|H（z）|在|z|>0处收敛，极点全部在z=0处，即FIR一定为稳定系统；结构上是非递归的，没有输出到输入的反馈。对于数字均衡器，一方面用于处理语音信号，另一方面需要用到频率特性分段的带通滤波器，因此本文选用IIR滤波器。3.2 IIR滤波器设计（1）设计模拟滤波器。根据实际需要确定滤波器的参数，利用的滤波器的设计公式设计模拟滤波器并得到其传递函数H（s）；（2）将滤波器转化为数字滤波器。利用冲激响应不变法或双线性变法将H（s）转化成H（z）；（3）频率变换。上述得到的是低通滤波器，利用变换公式作频率变换得到高通、带通滤波器。3.3均衡器的

18、原理均衡器的基本功能是调节各个频段内的信号的强弱，设计过程如下：（1） 设计出八个频段的对应的带通滤波器；（2） 把语音信号分8路用八个带通滤波器进行滤波；（3） 将八个滤波器的输出结果进行加权求和，yn=i=1Nkiyin *（25） ki表示第i路的权值，第i路的输出为 yi(n)，y(n)表示均衡器的输出。三 项目完成过程1. 语音识别以下为语音识别过程所用到的相关函数介绍，具体关系为：VQ.m调用mfcc.m获取单个录音的mel倒谱系数，接着mfcc.m调用melfb.m，完成语音信号从普通频谱到mel频谱的转换。在train.m文件调用VQ.m，即该文件用于训练码本，从而获得码本库；

19、在test.m文件调用distance.m计算训练录音与测试录音的mel倒谱系数的欧氏距离，如果最小距离小于一个阈值，则可判断该说话人是否在码本库内以及说话人是谁。1.1 VQ函数function r = VQ(d,k)e = .01;r = mean(d, 2);dpr = 10000;for i = 1:log2(k) r = r*(1+e), r*(1-e); while (1 = 1) z = distance(d, r); m,ind = min(z, , 2); t = 0; for j = 1:2i r(:, j) = mean(d(:, find(ind = j), 2); x

20、 = distance(d(:, find(ind = j), r(:, j); for q = 1:length(x) t = t + x(q); end end if (dpr - t)/t) < e) break; else dpr = t; end endend1.2 mfcc函数function r = mfcc(s, fs)m = 100;n = 256;l = length(s);nbFrame = floor(l - n) / m) + 1; for i = 1:n for j = 1:nbFrame M(i, j) = s(j - 1) * m) + i); enden

21、dh = hamming(n); /加窗函数M2 = diag(h) * M;for i = 1:nbFrame frame(:,i) = fft(M2(:, i); /逐帧进行快速傅立叶变换endt = n / 2;tmax = l / fs;m = melfb(20, n, fs); /通过Mel滤波器组得到Mel频谱 n2 = 1 + floor(n / 2); / 求每帧要进行变换的数据的行数z = m * abs(frame(1:n2, :).2; /求矩阵绝对值和的平方并与原始矩阵相乘。r = dct(log(z); / dct变换1.3 melfb函数function

22、m = melfb(p, n, fs)f0 = 700 / fs;fn2 = floor(n/2);lr = log(1 + 0.5/f0) / (p+1);bl = n * (f0 * (exp(0 1 p p+1 * lr) - 1);/直接转换为FFT的数字模型b1 = floor(bl(1) + 1;b2 = ceil(bl(2);b3 = floor(bl(3);b4 = min(fn2, ceil(bl(4) - 1;pf = log(1 + (b1:b4)/n/f0) / lr;fp = floor(pf);pm = pf - fp;r = fp(b2:b4) 1+fp(1:b3

23、);c = b2:b4 1:b3 + 1;v = 2 * 1-pm(b2:b4) pm(1:b3);m = sparse(r, c, v, p, 1+fn2);1.4 train函数function code = train(traindir, n)k = 16; for i = 1:n / 对数据库中的代码形成码本 file = sprintf('%svoice%d.wav', traindir, i); disp(file); s, fs = wavread(file); v = mfcc(s, fs); / 计算 MFCC's 提取特征特征，返回值是Mel倒谱系数

24、，是一个log的dct得到的 codei = VQ(v, k); / 训练VQ码本 通过矢量量化，得到原说话人的VQ码本end1.5 distance函数function d = distance(x, y)M, N = size(x);M2, P = size(y); /音频y赋值给M2，Pif (M = M2) error('不匹配！') /两个音频时间长度不相等endd = zeros(N, P);if (N < P) /在两个音频时间长度相等的前提下 copies = zeros(1,P); for n = 1:N d(n,:) = sum(x(:, n+copi

25、es) - y) .2, 1); endelse copies = zeros(1,N); for p = 1:P d(:,p) = sum(x - y(:, p+copies) .2, 1)' end /成对欧氏距离的两个矩阵的列之间的距离endd = d.0.5;1.6 test函数function finalmsg = test(testdir, n, code) nameList='胡程勇','周轶轩','刘杨雨','腾飞','夏文杰' file = sprintf('%stest.wav&

26、#39;, testdir);s, fs = wavread(file); v = mfcc(s,fs); / 计算MFCCdistmin = inf; /无穷k1 = 0;for i = 1:length(code) / code的最大长度，对于每个训练码书，计算失真d=distance(v,codei); dist = sum(min(d,2) / size(d,1); /size(d,1): 列的长度即行数 / min(A,dim),dim取1或2，取1时，与max(A)完全相同，取2时，该函数返回一个列向量， / min第i个元素是A矩阵的第i行上的最小值 if dist < d

27、istmin distmin = dist; k1=i; end enddisp(distmin);if(distmin < 4.5) xianshi = sprintf(' 说话者为 %s !', nameListk1);else xianshi = ' 说话者不在数据库中！ 'endmsgbox(xianshi);上述diatance.m函数中，阈值4.5是通过将码本库里的语音分别作为测试音，与码本库里的码本计算得到的最小距离，另外再用两个不是码本库里的测试音与码本库的码本计算得到最小距离，取这些最小距离的平均值作为此次识别的阈值，如果小于阈值，则可说

28、明该说话者是码本库中的一个，并进一步可以判断出是码本库中的哪一个。2. 噪声降除Wiener Scalart:NoiseCounter=0; % 初始化NoiseCounterNoiseLength=9; % 设置噪声平滑区间长度G=ones(size(N); % 初始化谱估计器Gamma=G;X=zeros(size(Y); % 初始化Xfor i=1:numberOfFrames SpeechFlag=SF(i); if i<=NIS % 若i<=NIS在前导无声(噪声)段 SpeechFlag=0; NoiseCounter=100; %else % i>NIS判断是否

29、为有话帧 %NoiseFlag, SpeechFlag, NoiseCounter, Dist=vad(Y(:,i),N,NoiseCounter); end if SpeechFlag=0 % 在无话段中平滑更新噪声谱值 N=(NoiseLength*N+Y(:,i)/(NoiseLength+1); LambdaD=(NoiseLength*LambdaD+(Y(:,i).2)./(1+NoiseLength);%更新和平滑噪声方差 end gammaNew=(Y(:,i).2)./LambdaD; % 计算后验信噪比 xi=alpha*(G.2).*Gamma+(1-alpha).*ma

30、x(gammaNew-1,0); % 计算先验信噪比 Gamma=gammaNew; G=(xi./(xi+1); % 计算维纳滤波器的谱估计器X(:,i)=G.*Y(:,i); % 维纳滤波后的幅值LMS adaptive filter:s=s/max(abs(s); % 幅值归一化N=length(s); % 求出信号长度global time;time=(0:N-1)/fs; % 设置时间ns=1.5*cos(2*pi*50*time); % 计算出50Hz工频信号global x;x=s+ns' % 语音信号和50Hz工频信号叠加snr1=SNR_singlech(s,x);

31、% 计算叠加50Hz工频信号后的信噪比x1=cos(2*pi*50*time); % 设置x1和x2x2=sin(2*pi*50*time);w1=0.1; % 初始化w1和w2w2=0.1;e=zeros(1, N); % 初始化e和yy=zeros(1, N);mu=0.05; % 设置mufor i=1: N % LMS自适应陷波器滤波 y(i)=w1 * x1(i)+ w2 * x2(i); e(i) =x(i)-y(i); w1=w1+mu * e(i) * x1(i); w2=w2+mu * e(i) * x2(i);endglobal outputoutput=e' %

32、陷波器输出snr2=SNR_singlech(s,output); % 计算滤波后的信噪比 snr=snr2-snr1;3. 语音均衡器该部分GUI界面如下，主要分成4个部分，音频选择及播放，音感特征调节以及均衡，均衡后音频播放，以及相关的波形图显示。3.1音频选择这里既可以选择音频文件，也可以现场录制一段语音进行均衡处理；3.2 音感特征调节及均衡音感特征调节实质是调节八个控件，从而改变对应滤波器的权值；滤波器类型选择，用户可以选择4种滤波器，分别为巴特沃斯滤波器，切比雪夫1滤波器，切比雪夫2滤波器以及椭圆滤波器。均衡器的功能的实现，其实质就是生成带通滤波器，并对原信号进行滤波。首先进行频段

33、划分handles.fband=20 100 200 500 1000 2000 4000 8000 16000;接着便是产生滤波器，该过程我们在确定选择按钮的回调函数中实现, for i= 1:8 if get(handles.popufilter,'Value') = 2 b a=butter(order,2*fband(i:i+1)/handles.Fs); /巴特沃斯滤波器（阶数，3db截止频率的归一化值） elseif get(handles.popufilter,'Value') = 3 b a=cheby1(order,0.5,2*fband(i:

34、i+1)/handles.Fs);/切比雪夫一型（阶数，通带最大衰减，通带截止频率归一化值） elseif get(handles.popufilter,'Value') = 4 b a=cheby2(order,20,2*fband(i:i+1)/handles.Fs);/切比雪夫二型（阶数，阻带最小衰减，通带截止频率归一化值） elseif get(handles.popufilter,'Value') = 5 b a=ellip(order,0.5,20,2*fband(i:i+1)/handles.Fs); /椭圆滤波器 end eval(sprintf

35、('k = get(handles.band%d,''Value'');',i); /获取相应的权值y =(9k-1)/8*filter(b,a,handles.y); /对语音信号进行滤波 handles.yy = handles.yy + y; /将滤波后的输出波形叠加 y =(9k-1)/8*filter(b,a,handles.y);的作用是加权以便下一句的求和，因为音感特征调节所获得的值在（0,2）之间，而需要的权值是（0,10）之间，所以令w为需要的权值，由于实际需要w为指数变化，故设代入w=0，k=0；w=1，k=1；w=10，k=

36、2，可以解得a=1/8，b=9，c=-1，从而得到上述的公式。3.3 均衡后播放对均衡过后的语音信号进行播放；3.4 相关波形的显示该部分主要显示均衡之前的语音信号的时域波形图，以及均衡之后的时域波形图，方便用户进行比较。4. GUI界面设计本项目通过MATLAB的GUI设计界面，使用户更加方便的运用这些功能，设计出来的主界面如下，分为三个功能，语音识别，降噪以及均衡器。语音识别界面如下：训练部分我们设计了用户可以自由选择码本训练的数量，也可以现场进行码本录制，然后进行码本训练；识别部分主要进行测试音的录制以及识别；界面的右部分主要是录音波形的显示。降噪界面如下： 我们分别对音频信号加入工频噪

37、声（这里的工频噪声是一个50Hz的余弦信号）和白噪声，并且针对这两种噪声使用不同的去噪方法进行解决，并且我们可以现场录音，将现场的录音中带有的噪声用维纳滤波器除去。可以看到我们给出了原始音频信号频谱，加噪信号频谱和去噪信号频谱，以此给读者一个更加直观的去噪过程。均衡器界面部分已在第三部分介绍，这里不再叙述。四 项目结果与分析1. 语音识别识别结果显示如下，该部分设计过程中，我们发现利用VQ矢量量化的方法存在一些问题，如果录音的时候不说话，那么识别判断出来的结果并不是显示“说话人不在语音库中”，而是上一次的判断结果，即如果上一次判断是库中的“胡程勇”，那么下一次录音时无人说话，则判断出来的还是库中的“胡程勇”。我们在多次实验中发现，一般将库中的语音作为测试音，所得到的最小距离都比较大，大概在4-5之间，但是如果是空白的录音，其测试所得到的最小