#matlab处理音频信号.docx

#matlab处理音频信号.docx

《#matlab处理音频信号.docx》由会员分享，可在线阅读，更多相关《#matlab处理音频信号.docx（16页珍藏版）》请在冰点文库上搜索。

#matlab处理音频信号

Matlab处理音频信号

一、问题的提出：

数字语音是信号的一种，我们处理数字语音信号，也就是对一种信号的处理，那信号是什么呢？

信号是传递信息的函数。

一、问题的提出：

数字语音是信号的一种，我们处理数字语音信号，也就是对一种信号的处理，那信号是什么呢？

信号是传递信息的函数。

离散时间信号%26mdash。

%26mdash。

序列%26mdash。

%26mdash。

可以用图形来表示。

按信号特点的不同，信号可表示成一个或几个独立变量的函数。

例如，图像信号就是空间位置<二元变量）的亮度函数。

一维变量可以是时间，也可以是其他参量，习惯上将其看成时间。

信号有以下几种：

<1）连续时间信号：

在连续时间范围内定义的信号，但信号的幅值可以是连续数值，也可以是离散数值。

当幅值为连续这一特点情况下又常称为模拟信号。

实际上连续时间信号与模拟信号常常通用，用以说明同一信号。

<2）离时间信号：

时间为离散变量的信号，即独立变量时间被量化了。

而幅度仍是连续变化的。

<3）数字信号：

时间离散而幅度量化的信号。

语音信号是基于时间轴上的一维数字信号，在这里主要是对语音信号进行频域上的分析。

在信号分析中，频域往往包含了更多的信息。

对于频域来说，大概有8种波形可以让我们分析：

矩形方波，锯齿波，梯形波，临界阻尼指数脉冲波形，三角波，余旋波，余旋平方波，高斯波。

对于各种波形，我们都可以用一种方法来分析，就是傅立叶变换：

将时域的波形转化到频域来分析。

于是，本课题就从频域的角度对信号进行分析，并通过分析频谱来设计出合适的滤波器。

当然，这些过程的实现都是在MATLAB软件上进行的，MATLAB软件在数字信号处理上发挥了相当大的优势。

二、设计方案：

利用MATLAB中的wavread命令来读入<采集）语音信号，将它赋值给某一向量。

再将该向量看作一个普通的信号，对其进行FFT变换实现频谱分析，再依据实际情况对它进行滤波。

对于波形图与频谱图<包括滤波前后的对比图）都可以用MATLAB画出。

我们还可以通过sound命令来对语音信号进行回放，以便在听觉上来感受声音的变化。

选择设计此方案，是对数字信号处理的一次实践。

在数字信号处理的课程学习过程中，我们过多的是理论学习，几乎没有进行实践方面的运用。

这个课题正好是对数字语音处理的一次有利实践，而且语音处理也可以说是信号处理在实际应用中很大众化的一方面。

这个方案用到的软件也是在数字信号处理中非常通用的一个软件%26mdash。

%26mdash。

MATLAB软件。

所以这个课题的设计过程也是一次数字信号处理在MATLAB中应用的学习过程。

课题用到了较多的MATLAB语句，而因为课题研究范围所限，真正与数字信号有关的命令函数却并不多。

三、主体部分：

<一）、语音的录入与打开：

[y,fs,bits]=wavread（'Blip',[N1N2]>。

用于读取语音，采样值放在向量y中，fs表示采样频率（Hz>，bits表示采样位数。

[N1N2]表示读取从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的FFT。

且和x相同长度。

若x为一矩阵，则y是对矩阵的每一列向量进行FFT。

如果x长度是2的幂次方，函数fft执行高速基－2FFT算法；否则fft执行一种混合基的离散傅立叶变换算法，计算速度较慢。

函数FFT的另一种调用格式为y=fft（x,N>

式中，x，y意义同前，N为正整数。

函数执行N点的FFT。

若x为向量且长度小于N，则函数将x补零至长度N。

若向量x的长度大于N，则函数截短x使之长度为N。

若x为矩阵，按相同方法对x进行处理。

经函数fft求得的序列y一般是复序列，通常要求其幅值和相位。

MATLAB提供求复数的幅值和相位函数：

abs，angle，这些函数一般和FFT同时使用。

函数abs（x>用于计算复向量x的幅值，函数angle（x>用于计算复向量的相角，介于和之间，以弧度表示。

函数unwrap（p>用于展开弧度相位角p,当相位角绝对变化超过时，函数把它扩展至。

用MATLAB工具箱函数fft进行频谱分析时需注意：

<1）函数fft返回值y的数据结构对称性

若已知序列x=[4,3,2,6,7,8,9,0],求X（k>=DFT[x（n>]。

利用函数fft计算，用MATLAB编程如下：

N=8。

n=0:

N-1。

xn=[43267890]'。

XK=fft（xn>

结果为：

XK=

39.0000

-10.7782+6.2929i

0-5.0000i

4.7782-7.7071i

5.0000

4.7782+7.7071i

0+5.0000i

-10.7782-6.2929i

由程序运行所得结果可见，X（k>和x（n>的维数相同，共有8个元素。

X（k>的第一行元素对应频率值为0，第五行元素对应频率值为Nyquist频率，即标准频率为1.因此第一行至第五行对应的标准频率为0～1。

而第五行至第八行对应的是负频率，其X（k>值是以Nyquist频率为轴对称。

<注：

通常表示为Nyquist频率外扩展，标以正值。

）

一般而言，对于N点的x（n>序列的FFT是N点的复数序列，其点n=N/2+1对应Nyquist频率，作频谱分析时仅取序列X（k>的前一半，即前N/2点即可。

X（k>的后一半序列和前一半序列时对称的。

<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',[10245120]>。

sound（x,fs,bits>。

X=fft（x,4096>。

magX=abs（X>。

angX=angle（X>。

subplot（221>。

plot（x>。

title（'原始信号波形'>。

subplot（222>。

plot（X>。

title（'原始信号频谱'>。

subplot（223>。

plot（magX>。

title（'原始信号幅值'>。

subplot（224>。

plot（angX>。

title（'原始信号相位'>。

程序运行可以听到声音，得到的图形为：

<2）定点分析：

已知一个语音信号，数据采样频率为100Hz，试分别绘制N＝128点DFT的幅频图和N＝1024点DFT幅频图。

编程如下：

x=wavread（'ding.wav'>。

sound（x>。

fs=100。

N=128。

y=fft（x,N>。

magy=abs（y>。

f=（0:

length（y>-1>'*fs/length（y>。

subplot（221>。

plot（f,magy>。

xlabel（'频率。

ylabel（'幅值'>。

title（'N=128（a>'>。

grid

subplot（222>。

plot（f（1:

N/2>,magy（1:

N/2>>。

xlabel（'频率。

ylabel（'幅值'>。

title（'N=128（b>'>。

grid

fs=100。

N=1024。

y=fft（x,N>。

magy=abs（y>。

f=（0:

length（y>-1>'*fs/length（y>。

subplot（223>。

plot（f,magy>。

xlabel（'频率。

ylabel（'幅值'>。

title（'N=1024（c>'>。

grid

subplot（224>。

plot（f（1:

N/2>,magy（1:

N/2>>。

xlabel（'频率。

ylabel（'幅值'>。

title（'N=1024（d>'>。

grid

运行结果如图：

上图（a>、（b>为N=128点幅频谱图，（c>、（d>为N=1024点幅频谱图。

因为采样频率f=100Hz，故Nyquist频率为50Hz。

（a>、（c>是0～100Hz频谱图，（b>、（d>是0～50Hz频谱图。

由（a>或（c>可见，整个频谱图是以Nyquist频率为轴对称的。

因此利用fft对信号作频谱分析，只要考察0～Nyquist频率<采样频率一半）范围的幅频特性。

比较（a>和（c>或（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>>。

grid

xlabel（'Hz'>,ylabel（'|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>>。

grid

xlabel（'Hz'>,ylabel（'|H（ejw>|'>

结果如图：

<三）、滤波器设计：

1、相关原理：

设计数字滤波器的任务就是寻求一个因果稳定的线性时不变系统，并使系统函数H（z>具有指定的频率特性。

数字滤波器从实现的网络结构或者从单位冲激响应分类，可以分成无限长单位冲激响应

数字滤波器频率响应的三个参数：

<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）卷积函数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中的数据进行滤波，结果放入向量y。

<3）函数fftfilt

函数fftfilt的调用格式为y=fftfilt（b,x>

该格式是利用基于FFT的重叠相加法对数据进行滤波，这种频域滤波技术只对FIR滤波器有效。

该函数是通过向量b描述的滤波器对x数据进行滤波。

关于用butter函数求系统函数分子与分母系数的几种形式。

[b,a]=butter（N,wc,'high'>:

设计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%26rdquo。

，表示设计的是模拟滤波器。

这里不作讨论。

为了设计任意的选项巴特沃斯滤波器，必须知道阶数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。

ws1%26lt。

ws2%26lt。

wp2。

2、设计内容：

<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=filter（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。

[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（'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>。

title（'IIR滤波后信号的频谱'>。

z=fftfilt（b,x>。

Z=fft（z>。

subplot（325>。

plot（z>。

title（'FIR滤波后信号的波形'>。

subplot（326>。

plot（Z>。

title（'FIR滤波后信号的频谱'>。

得到结果如图：

<4）2N阶带通滤波器的设计<在这里，以10阶为例，其中wc为其3dB边缘频率，以为单位，wc=[w1,w2],w1wcw2），程序设计如下：

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（'滤波前信号的波形'>。

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（'FIR滤波后信号的频谱'>。

得到结果如图：

<5）2N阶带阻滤波器的设计<在这里，以10阶为例，其中wc为其3dB边缘频率，以为单位，wc=[w1,w2],w1wcw2），程序设计如下：

x=wavread（'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滤波后信号的波形'>。

subplot（326>。

plot（Z>。

title（'FIR滤波后信号的频谱'>。

得到结果如图：

<6）小结：

以上几种滤波，我们都可以从信号滤波前后的波形图以及频谱图上看出变化。

当然，也可以用sound（>函数来播放滤波后的语音，从听觉上直接感受语音信号的变化，但因为人耳听力的限制，有些情况下我们是很难听出异同的。

同样，通过函数的调用，也可以将信号的频谱进行%26ldquo。

分离观察%26rdquo。

，如显出信号的幅值或相位。

下面，通过改变系统函数的分子与分母系数比，来观察信号滤波前后的幅值与相位。

并且使结果更加明显，使人耳得以很容易的辨听。

x=wavread（'ding.wav'>。

sound（x>。

b=100。

a=5。

y=filter（b,a,x>。

X=fft（x,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>>。

title（'滤波前信号相位'>。

subplot（212>。

plot（angle（Y>>。

title（'滤波后信号相位'>。

得图：

可以看到相位谱没什么变化。

<四）、界面设计：

直接用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.为每个回调函数添加代码以实现GUI程序的具体功能。

这一步与一般函数文件的编辑调试过程相同。

设计过程及内容：

在MATLAB版面上，通过键入GUIDE弹出一个菜单栏进入gui制作界面<或者在File到new来进入gui）,从而开始应用界面的制作。

该界面主要实现了以下几个功能：

①打开wav格式的音频文件，并将该音频信号的值读取并赋予某一向量；

②播放音频文件，可以选择性的显示该音频信号的波形、频谱、幅值以及相位；

③对音频信号进行IIR与FIR的5阶固定滤波处理，可以选择性的显示滤波前后信号的波形、频谱、幅值以及相位，以及播放滤波后的声音。

界面如图所示：

通过该界面，可以方便用户进行语音信号的处理。

界面主程序见附件。

<五）、校验：

1、本设计圆满的完成了对语音信号的读取与打开，与课题的要求十分相符；

2、本设计也较好的完成了对语音信号的频谱分析，通过fft变换，得出了语音信号的频谱图；

3、在滤波这一块，课题主要是从巴特沃斯滤波器入手来设计滤波器，也从一方面基本实现了滤波；

4、初略的完成了界面的设计，但也存在相当的不足，只是很勉强的达到了打开语音文件、显示已定滤波前后的波形等图。

四、结论：

语音信号处理是语音学与数字信号处理技术相结合的交叉学科，课题在这里不讨论语音学，而