如何利用matlab处理音频信号

1、Matlab处理音频信号一、问题的提出：数字语音是信号的一种，我们处理数字语音信号，也就是对一 种信号的处理，那信号是什么呢？ 信号是传递信息的函数。一、问题的提出：数字语音是信号的一种，我们处理数字语音信号，也就是对一种信号的处理，那信 号是什么呢？信号是传递信息的函数。离散时间信号 26mdash;%26mdas市 列26mdash;%26mdash；以用图形来表示。按信号特点的不同，信号可表示成一个或几个独立变量的函数。例如，图像信号就是空间位置(二元变量)的亮度函数。一维变量可以是时间，也可以是其他参量，习惯上将其看成时间。信号有以下几种：(1)连续时间信号：在连续时间范围内定义的信号

2、，但信号的幅值可以是连续数值，也可以是离散数值。当幅值为连续这一特点情况下又常称为模拟信号。实际上 连续时间信号与模拟信号常常通用，用以说明同一信号。(2)离时间信号：时间为离散变量的信号，即独立变量时间被量化了。而幅度仍 是连续变化的。(3)数字信号：时间离散而幅度量化的信号。语音信号是基于时间轴上的一维数字信号，在这里主要是对语音信号进行频域上的分析。在信号分析中，频域往往包含了更多的信息。对于频域来说，大概有8种波形可以让我们分析：矩形方波，锯齿波，梯形波，临界阻尼指数脉冲波形，三角 波，余旋波，余旋平方波，高斯波。对于各种波形，我们都可以用一种方法来分 析，就是傅立叶变换：将时域的波形

3、转化到频域来分析。于是，本课题就从频域的角度对信号进行分析，并通过分析频谱来设计出合适的滤 波器。当然,这些过程的实现都是在 MATLAB软件上进行的，MATLAB软件在数 字信号处理上发挥了相当大的优势。二、设计方案：利用MATLAB中的wavrea而令来读入（采集）语音信号，将它赋值给某一向量。再将该向量看作一个普通的信号，对其进行 FFT变换实现频谱分析，再依据 实际情况对它进行滤波。对于波形图与频谱图（包括滤波前后的对比图）都可以用 MATLAB画出。我们还可以通过sound命令来对语音信号进行回放，以便在听觉 上来感受声音的变化。选择设计此方案，是对数字信号处理的一次实践。在数字信号

4、处理的课程学习过程 中，我们过多的是理论学习，几乎没有进行实践方面的运用。这个课题正好是对数 字语音处理的一次有利实践，而且语音处理也可以说是信号处理在实际应用中很大 众化的一方面。这个方案用到的软件也是在数字信号处理中非常通用的一个软件26mdash;%26mdash;MATLA啾件。所以这个课题的设计过程也是一次数字信号处理在MATLAB中应用的学习过程。课题用到了较多的 MATLAB语句，而由于 课题研究范围所限，真正与数字信号有关的命令函数却并不多。三、主体部分：（一）、语音的录入与打开：y,fs,bits=wavread（Blip,N1 N2）用于读取语音，采样值放在向量 y中，fs

5、表示采样 频率（Hz）, bits表示采样位数。N1 N2表示读取从N1点到N2点的值（若只有一 个N的点则表示读取前N点的采样值）。sound（x,fs,bits）用于对声音的回放。向量y则就代表了一个信号（也即一个复杂 的26ldquo函数表达式26rdquo）也就是说可以像处理一个信号表达式一样处理 这个声音信号。FFT的MATLAB实现在MATLAB的信号处理工具箱中函数FFT和IFFT用于快速傅立叶变换和逆变 换。下面介绍这些函数。函数FFT用于序列快速傅立叶变换 函数的一种调用格式为y=fft(x)其中，x是序列，y是序列的FFT, x可以为一向量或矩阵，若x为一向量，y是x 的F

6、FT。且和x相同长度。若x为一矩阵，则y是对矩阵的每一列向量进行FFT。如果x长度是2的幕次方，函数fft执行高速基2FFT算法；否则fft执行一种混 合基的离散傅立叶变换算法，计算速度较慢。函数FFT的另一种调用格式为y=fft(x,N)式中，x, y意义同前，N为正整数。函数执行N点的FFT。若x为向量且长度小于N,则函数将x补零至长度No若 向量x的长度大于N,则函数截短x使之长度为No若x为矩阵，按相同方法对x 进行处理。经函数fft求得的序列y 一般是复序列，通常要求其幅值和相位。MATLAB提供求复数的幅值和相位函数：abs angle这些函数一般和FFT同时使用。函数abs(x用

7、于计算复向量x的幅值，函数angle(x用于计算复向量的相角，介于 和之间，以弧度表示。函数unwrap(p归于展开弧度相位角p，当相位角绝对变化超过时，函数把它扩展 至。用MATLAB工具箱函数fft进行频谱分析时需注意：(1)函数fft返回值y的数据结构对称性若已知序列 x=4,3,2,6,7,8,9,0胀 X(k尸DFTx(n)。利用函数fft计算，用MATLAB编程如下：N=8;n=0:N-1;xn=4 3 2 6 7 8 9 0;XK=fft(xn)结果为：XK =39.0000-10.7782 + 6.2929i0 - 5.0000i4.7782 - 7.7071i5.00004.

8、7782 + 7.7071i0 + 5.0000i-10.7782 - 6.2929i由程序运行所得结果可见，X(k)和x(n)的维数相同，共有8个元素。X(k)的第一行 元素对应频率值为0,第五行元素对应频率值为 Nyquist频率，即标准频率为1.因 此第一行至第五行对应的标准频率为01。而第五行至第八行对应的是负频率，其X(k)值是以Nyquist频率为轴对称。(注：通常表示为Nyquist频率外扩展，标以正值。)一般而言，对于N点的x(n)序列的FFT是N点的复数序列，其点n=N/2+1对应 Nyquist频率，作频谱分析时仅取序列 X(k)的前一半，即前N/2点即可。X(k)的后 一

9、半序列和前一半序列时对称的。(2)频率计算若N点序列x(n)(n=0,1，,N-1)是在采样频率 下获得的。它的FFT也是N点序列, 即X(k)(k=0,1,2,N-1),则第k点所对应实际频率值为f=k*f/N.(3)作FFT分析时，幅值大小与FFT选择点数有关，但不影响分析结果。2、设计内容：(1)下面的一段程序是语音信号在MATLAB中的最简单表现，它实现了语音的读入打开，以及绘出了语音信号的波形频谱图。x,fs,bits=wavread(ding.wav,1024 5120);sound(x,fs,bits);X=fft(x,4096);magX=abs(X);angX=angle(X

10、);subplot(221);plot(x);title媒始信号波形);subplot(222);plot(X); title(原始信号频谱);subplot(223);plot(magX);title原始信号幅值);subplot(224);plot(angX);title原始信号相位);程序运行可以听到声音，得到的图形为：100Hz,试分别绘制 N = 128(2)定点分析：已知一个语音信号，数据采样频率为 点DFT的幅频图和N = 1024点DFT幅频图。编程如下:x=wavread(ding.wav);sound(x);fs=100;N=128;y=fft(x,N);magy=abs(

11、y);f=(0:length(y)-1)*fs/length(y);subplot(221);plot(f,magy);xlabel(频率(Hz) );ylabel(fi值)；title(N=128(a);gridsubplot(222);plot(f(1:N/2),magy(1:N/2);xlabel(频率(Hz) );ylabel(fi值)；title(N=128(b);gridfs=100;N=1024;y=fft(x,N);magy=abs(y);f=(0:length(y)-1)*fs/length(y);subplot(223);plot(f,magy);xlabel（频率（Hz）

12、 ）;ylabel（|fi值）;title(N=1024(c);gridsubplot(224);plot(f(1:N/2),magy(1:N/2);xlabel(频率(Hz) );ylabel(lfi值)；title(N=1024(d);grid运行结果如图：上图(a) (b)为N=128点幅频谱图，(c)、(d)为N=1024点幅频谱图。由于采样频率f =100Hz,故 Nyquist频率为 50Hz。(a)、(c)是 0100Hz频谱图，(b)、(d)是 050Hz 频谱图。由(a域(c)可见，整个频谱图是以Nyquist频率为轴对称的。因此利用fft对 信号作频谱分析，只要考察0Nyq

13、uist频率(采样频率一半)范围的幅频特性。比较(a和或(b)和(d)可见，幅值大小与fft选用点数N有关，但只要点数N足够不影响研究结果。从上图幅频谱可见，信号中包括15Hz和40Hz的正弦分量。(3)若信号长度T=25.6s,即抽样后x(n)点数为T/Ts=256 ,所得频率分辨率为Hz,以此观察数据长度N的变化对DTFT分辨率的影响：编程如下：x,fs,bits=wavread(ding.wav);N=256;f=0:fs/N:fs/2-1/N;X=fft(x);X=abs(X);subplot(211)plot(f(45:60),X(45:60);gridxlabel(Hz),ylab

14、el(|H(ejw)|)%数据长度N扩大4倍后观察信号频谱N=N*4;f=0:fs/N:fs/2-1/N;X=fft(x);X=abs(X);subplot(212)plot(f(45*4:4*60),X(4*45:4*60);gridxlabel(Hz),ylabel(|H(ejw)|)结果如图：(三)、滤波器设计：1、相关原理：设计数字滤波器的任务就是寻求一个因果稳定的线性时不变系统，并使系统函数 H(z)具有指定的频率特性。数字滤波器从实现的网络结构或者从单位冲激响应分类，可以分成无限长单位冲激 响应(IIR)数字滤波器和有限长单位冲激响应(FIR)数字滤波器。数字滤波器频率响应的三个参

15、数：(1)幅度平方响应:(2)相位响应其中，相位响应(3)群时延响应IIR数字滤波器：IIR数字滤波器的系统函数为的有理分数，即IIR数字滤波器的逼近问题就是求解滤波器的系数 和，使得在规定的物理意义上 逼近所要求的特性的问题。如果是在 s平面上逼近，就得到模拟滤波器，如果是在 z平面上逼近，则得到数字滤波器。FIR数字滤波器：设FIR的单位脉冲响应h(n)为实数，长度为N,则其z变换和频率响应分别为按频域采样定理FIR数字滤波器的传输函数 H(z)和单位脉冲响应h(n)可由它的N 个频域采样值H(k)唯一确定。MATLAB中提供了几个函数，分别用于实现IIR滤波器和FIR滤波器。(1)卷积函

16、数conv卷积函数conv的调用格式为c=conv(a,b)该格式可以计算两向量a和b的卷积，可以直接用于对有限长信号采用FIR滤波器的滤波。(2)函数 filter函数filter的调用格式为y=filter(b,a,x)该格式采用数字滤波器对数据进行滤波，既可以用于 IIR滤波器，也可以用于FIR 滤波器。其中向量b和a分别表示系统函数的分子、分母多项式的系数，若a=1,此时表示FIR滤波器，否则就是IIR滤波器。该函数是利用给出的向量b和a,对x中的数据进行滤波，结果放入向量 V。3 3)函数 fftfilt函数fftfilt的调用格式为y=fftfilt(b,x)该格式是利用基于FFT

17、的重叠相加法对数据进行滤波，这种频域滤波技术只对 FIR滤波器有效。该函数是通过向量 b描述的滤波器对x数据进行滤波。关于用butter函数求系统函数分子与分母系数的几种形式。b,a=butter(N,wGhigh)设计N阶高通滤波器，wc为它的3dB边缘频率，以 为单 位，故。b,a=butter(N,wc):当wc为具有两个元素的矢量 wc=w1,w2时，它设计2N阶带通滤 波器，3dB通带为，w的单位为。b,a=butter(N,wc,stop)若 wc=w1,w2,则它设计 2N 阶带阻滤波器，3dB通带为， w的单位为。如果在这个函数输入变元的最后，加一个变元26ldquo;s%26

18、rdquq;表示设计的是模拟滤波器。这里不作讨论。为了设计任意的选项巴特沃斯滤波器，必须知道阶数 N和3dB边缘频率矢量wc。 这可以直接利用信号处理工具箱中的 buttord函数来计算。如果已知滤波器指 标，和，则调用格式为N,wc=buttord(wp,ws,Rp,As)对于不同类型的滤波器，参数 wp和ws有一些限制：对于低通滤波器，wp%26lt;ws；对于高通滤波器，wp%26gt;ws；对于带通滤波器，wp和ws分别为具有两个元素的矢量，wp=wp1,wp2和ws=ws1,ws2,并且ws1%26lt;wp1%26lt;wp2%26lt;ws2 对于带阻滤波器wp1%26lt;ws

19、1%26lt;ws2%26lt;wp22、设计内容:(1)滤波器示例：在这里为了说明如何用 MATLAB来实现滤波，特举出一个简单的函数信号滤波实 例(对信号x(n尸sin( n/4)+5cos( n/2)进行滤波，信号长度为500点)，从中了解滤 波的实现过程。程序如下：Wn=0.2*pi;N=5;b,a=butter(N,Wn/pi);n=0:499;x=sin(pi*n/4)+5*cos(pi*n/2);X=fft(x,4096);subplot(221);plot(x);title波前信号的波形);subplot(222);plot(X);title(据波前信号的频谱);y=filte

20、r(b,a,x);Y=fft(y,4096);subplot(223);plot(y);title就波后信号的波形);subplot(224);plot(Y);title波后信号的频谱);在这里，是采用了 butter命令，设计出一个巴特沃斯低通滤波器，从频谱图中可以 很明显的看出来。下面，也就是本课题的主要内容，也都是运用到了butter函数,以便容易的得到系统函数的分子与分母系数，最终以此来实现信号的滤波。(2) N阶高通滤波器的设计(在这里，以 5阶为例，其中wc为其3dB边缘频 率，以为单位)，程序设计如下：x=wavread(ding.wav);sound(x);N=5;wc=0.3

21、;b,a=butter(N,wc,high);X=fft(x);subplot(321);plot(x);title波前信号的波形);subplot(322);plot(X);title(据波前信号的频谱);y=filter(b,a,x);Y=fft(y);subplot(323);plot(y);title(IIR滤波后信号的波形);subplot(324);plot(Y);title(IIR滤波后信号的频谱);z=fftfilt(b,x);Z=fft(z);subplot(325);plot(z);title(FIR滤波后信号的波形);subplot(326);plot(Z);title(

22、FIR滤波后信号的频谱);得到结果如图:(3) N阶低通滤波器的设计(在这里，同样以5阶为例，其中wc为其3dB边缘频率，以为单位)，程序设计如下：x=wavread(ding.wav);sound(x);N=5;wc=0.3;b,a=butter(N,wc);X=fft(x);subplot(321);plot(x);title波前信号的波形);subplot(322);plot(X);title(据波前信号的频谱);y=filter(b,a,x);Y=fft(y);subplot(323);plot(y);title(IIR滤波后信号的波形);subplot(324);plot(Y);ti

23、tle(IIR滤波后信号的频谱);z=fffilt(b,x);Z=fft(z);subplot(325);plot(z);title(FIR滤波后信号的波形);subplot(325);plot(z);title(FIR滤波后信号的波形);得到结果如图：(4) 2N阶带通滤波器的设计(在这里，以10阶为例，其中wc为其3dB边缘频 率，以 为单位，wc=w1,w2,w1 wc w2 ,程序设计如下：x=wavread(ding.wav);sound(x);N=5;wc=0.3,0.6;b,a=butter(N,wc);X=fft(x);subplot(321);plot(x);title波前信

24、号的波形);subplot(322);plot(X);title(据波前信号的频谱);y=filter(b,a,x);Y=fft(y);subplot(323);plot(y);title(IIR滤波后信号的波形);subplot(324);plot(Y);title(IIR滤波后信号的频谱);z=ff情lt(b,x);Z=fft(z);subplot(326);plot(Z);title(FIR滤波后信号的频谱);得到结果如图：(5) 2N阶带阻滤波器的设计(在这里，以10阶为例，其中wc为其3dB边缘频率，以 为单位，wc=w1,w2,w1 wc w2 ,程序设计如下：x=wavread(

25、ding.wav);sound(x);N=5;wc=0.2,0.7;b,a=butter(N,wc,stop);X=fft(x);subplot(321);plot(x);title波前信号的波形);subplot(322);plot(X);title(据波前信号的频谱);y=filter(b,a,x);Y=fft(y);subplot(323);plot(y);title(IIR滤波后信号的波形);subplot(324);plot(Y);title(IIR滤波后信号的频谱);z=fftfilt(b,x);Z=fft(z);subplot(325);plot(z);title(FIR滤波后信

26、号的波形);得到结果如图：(6)小结：以上几种滤波，我们都可以从信号滤波前后的波形图以及频谱图上看 出变化。当然，也可以用sounder数来播放滤波后的语音，从听觉上直接感受语 音信号的变化，但由于人耳听力的限制，有些情况下我们是很难听出异同的。同样，通过函数的调用，也可以将信号的频谱进行26ldquo分离观察26rdquo;,如显出信号的幅值或相位。下面，通过改变系统函数的分子与分母系数比，来观察 信号滤波前后的幅值与相位。并且使结果更加明显，使人耳得以很容易的辨听。x=wavread(ding.wav);sound(x);b=100;a=5;y=filter(b,a,x);X=fft(x,

27、4096);subplot(221);plot(x);title波前信号的波形);subplot(222);plot(abs(X);title谍波前信号的幅值);Y=fft(y,4096);subplot(223);plot(y);title(*波后信号的波形);subplot(224);plot(abs(Y);title滤波后信号的幅值);结果如图：%26gt;%26gt; sound(y);可以听到声音明显变得高亢了。从上面的波形与幅值(即幅频)图，也可看出，滤 波后的幅值变成了滤波前的20倍。%26gt;%26gt; figure,subplot(211);plot(angle(X);t

28、itle滤波前信号相位);subplot(212);plot(angle(Y);title(g波后信号相位);得图：可以看到相位谱没什么变化。(四)、界面设计：直接用M文件编写GUI程序很繁琐，而使用GUIDE设计工具可以大大提高工作 效率。GUIDE相当于一个控制面板，从中可以调用各种设计工具以辅助完成界面 设计任务，例如控件的创建和布局、控件属性的编辑和菜单设计等。使用GUIDE设计GUI程序的一般步骤如下：1 .将所需控件从控件面板拖拽到 GUIDE的设计区域；2 .利用工具条中的工具(或相应的菜单和现场菜单)，快速完成界面布局；3 .设置控件的属性。尤其是tag属性，它是控件在程序内部的唯一标识；4 .如果需要，打开菜单编辑器为界面添加菜单或现场菜单；5 .保存设计。GUIDE默认把GUI程序保存为两个同名文件：一个是.fig文件，用 来保存窗体布局和所有控件的界面信息；一个是.m文件，该文件的初始内容是 GUIDE自动产生的程序框架，其中包括了各个控件回调函数的定义。该 M文件与一般的M文件没有本质区别，但是鉴于它的特殊性，MATALAB把这类文件统称为GUI-M文件。保存完后GUI-M文件自动在编辑调试器中打开以供编辑。6 .为每个回调函数添加代码以实现 G