matlab音频处理.docx
- 文档编号:27315905
- 上传时间:2023-06-29
- 格式:DOCX
- 页数:16
- 大小:22.64KB
matlab音频处理.docx
《matlab音频处理.docx》由会员分享,可在线阅读,更多相关《matlab音频处理.docx(16页珍藏版)》请在冰豆网上搜索。
matlab音频处理
matlab处理音频信号
一、 问题的提出:
数字语音是信号的一种,我们处理数字语音信号,也就是对一种信号的处理,那信号是什么呢?
信号是传递信息的函数.离散时间信号——序列—-可以用图形来表示。
按信号特点的不同,信号可表示成一个或几个独立变量的函数。
例如,图像信号就是空间位置(二元变量)的亮度函数。
一维变量可以是时间,也可以是其他参量,习惯上将其看成时间。
信号有以下几种:
(1)连续时间信号:
在连续时间范围内定义的信号,但信号的幅值可以是连续数值,也可以是离散数值.当幅值为连续这一特点情况下又常称为模拟信号.实际上连续时间信号与模拟信号常常通用,用以说明同一信号.
(2)离时间信号:
时间为离散变量的信号,即独立变量时间被量化了.而幅度仍是连续变化的。
(3)数字信号:
时间离散而幅度量化的信号。
语音信号是基于时间轴上的一维数字信号,在这里主要是对语音信号进行频域上的分析。
在信号分析中,频域往往包含了更多的信息。
对于频域来说,大概有8种波形可以让我们分析:
矩形方波,锯齿波,梯形波,临界阻尼指数脉冲波形,三角波,余旋波,余旋平方波,高斯波。
对于各种波形,我们都可以用一种方法来分析,就是傅立叶变换:
将时域的波形转化到频域来分析。
于是,本课题就从频域的角度对信号进行分析,并通过分析频谱来设计出合适的滤波器。
当然,这些过程的实现都是在MATLAB软件上进行的,MATLAB软件在数字信号处理上发挥了相当大的优势。
二、 设计方案:
利用MATLAB中的wavread命令来读入(采集)语音信号,将它赋值给某一向量.再将该向量看作一个普通的信号,对其进行FFT变换实现频谱分析,再依据实际情况对它进行滤波。
对于波形图与频谱图(包括滤波前后的对比图)都可以用MATLAB画出。
我们还可以通过sound命令来对语音信号进行回放,以便在听觉上来感受声音的变化。
选择设计此方案,是对数字信号处理的一次实践.在数字信号处理的课程学习过程中,我们过多的是理论学习,几乎没有进行实践方面的运用。
这个课题正好是对数字语音处理的一次有利实践,而且语音处理也可以说是信号处理在实际应用中很大众化的一方面.
这个方案用到的软件也是在数字信号处理中非常通用的一个软件——MATLAB软件。
所以这个课题的设计过程也是一次数字信号处理在MATLAB中应用的学习过程。
课题用到了较多的MATLAB语句,而由于课题研究范围所限,真正与数字信号有关的命令函数却并不多。
三、 主体部分:
(一)、语音的录入与打开:
[y,fs,bits]=wavread(’Blip’,[N1N2]);用于读取语音,采样值放在向量y中,fs表示采样频率(Hz),bits表示采样位数。
[N1N2]表示读取从N1点到N2点的值(若只有一个N的点则表示读取前N点的采样值)。
sound(x,fs,bits);用于对声音的回放。
向量y则就代表了一个信号(也即一个复杂的“函数表达式”)也就是说可以像处理一个信号表达式一样处理这个声音信号。
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('频率(Hz)');ylabel(’幅值’);
title(’N=128(a)’);grid
subplot(222);plot(f(1:
N/2),magy(1:
N/2));
xlabel(’频率(Hz)');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(’频率(Hz)');ylabel(’幅值’);
title('N=1024(c)');grid
subplot(224);plot(f(1:
N/2),magy(1:
N/2));
xlabel('频率(Hz)’);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)具有指定的频率特性。
数字滤波器从实现的网络结构或者从单位冲激响应分类,可以分成无限长单位冲激响应(IIR)数字滤波器和有限长单位冲激响应(FIR)数字滤波器.
数字滤波器频率响应的三个参数:
(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的单位为。
如果在这个函数输入变元的最后,加一个变元“s”,表示设计的是模拟滤波器。
这里不作讨论.
为了设计任意的选项巴特沃斯滤波器,必须知道阶数N和3dB边缘频率矢量wc。
这可以直接利用信号处理工具箱中的buttord函数来计算.如果已知滤波器指标,,和,则调用格式为
[N,wc]=buttord(wp,ws,Rp,As)
对于不同类型的滤波器,参数wp和ws有一些限制:
对于低通滤波器,wp 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()函数来播放滤波后的语音,从听觉上直接感受语音信号的变化,但由于人耳听力的限制,有些情况下我们是很难听出异同的. 同样,通过函数的调用,也可以将信号的频谱进行“分离观察”,如显出信号的幅值或相位.下面,通过改变系统函数的分子与分母系数比,来观察信号滤波前后的幅值与相位.并且使结果更加明显,使人耳得以很容易的辨听. 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(’滤波后信号的幅值'); 结果如图: 〉〉sound(y); 可以听到声音明显变得高亢了。 从上面的波形与幅值(即幅频)图,也可看出,滤波后的幅值变成了滤波前的20倍。 >>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、初略的完成了界面的设计,但也存在相当的不足,只是很勉强的达到了打开语音文件、显示已定滤波前后的波形等图。 四、结论: 语音信号处理是语音学与数字信号处理技术相结合的交叉学科,课题在这里不讨论语音学,而是将语音当做一种特殊的信号,即一种
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- matlab 音频 处理
![提示](https://static.bdocx.com/images/bang_tan.gif)