音频信号的采样与重构等matlab代码数字信号处理.docx
《音频信号的采样与重构等matlab代码数字信号处理.docx》由会员分享，可在线阅读，更多相关《音频信号的采样与重构等matlab代码数字信号处理.docx（17页珍藏版）》请在冰豆网上搜索。
音频信号的采样与重构等matlab代码数字信号处理
基于MATLAB的音频信号分析与处理
摘要
数字信号处理是一门发展十分迅速、应用非常广泛的前沿学科。
它的理论性和实践性都非常强。
MATLAB强大的计算仿真功能在数字信号处理领域起着非常大的作用。
出于对数字音频处理的兴趣，本文中将尝试利用所学的知识，如采样、滤波、重构等知识，对语音信号或是音频信号进行一定的处理。
本文详细给出了利用MATLAB对音频信号进行谱分析，信号滤波和重构的过程，加深了对所学数字信号处理知识的了解。
关键词：
滤波重构谱分析
Abstract
Digitalsignalprocessingisanadvancedsubjectwhichisquicklydevelopingandwidelyused.Itlaysagreatemphasisbothintheoryandpractice.MATLAB,thepowerfulcomputationandsimulationsoftware,playsagreatroleindigitalsignalprocessingfield.Asfortheinterestforthedigitalaudioprocessing,throughthispaper,Iamtryingtomakesomeprocessingaboutthesound（audio）signalwithwhatIhavelearntinclassroom,suchassampling,filtering,reconstruction,andsoon.Inthispaper,theprocessesoffrequencyamplitudeanalysis,filteringandconstructionoftheaudiosignal,whicharebasedonMATLABaredetailedlypresented,andIgainmoreunderstatingaboutknowledgeofdigitaldataprocessing.
Keywords:
filtering,reconstruction,frequencyamplitudeanalysis
摘要
Abstract
1数字滤波器
1.1数字滤波器概述
1.2IIR数字滤波器的设计理论
1.3用窗函数设计FIR滤波器
1.3.1设计思想
1.3.2．典型的窗函数
1.3.3．设计步骤
2快速傅立叶变换（FFT）
2.1FFT算法
2.2FFT的优越性
2.3用FFT进行频谱分析
3基于MATLAB的语音信号分析和处理
3.1MATLAB简介
3.2基于MATLAB的语音信号分析
3.2.1语音信号的采集及采样
3.2.2语音信号的频谱分析
3.2.3设计滤波器进行滤波处理
3.3基于MATLAB的语音信号的谱分析和重构
总结
参考文献
1数字滤波器
1.1数字滤波器概述
数字滤波器是对数字信号实现滤波的线性时不变系统。
它通过一定的计算或判断程序来减少干扰信号的比重，它实现的是程序滤波。
数字滤波器的输入与输出均为数字信号，通过一定的运算关系改变输入信号所含频率成分的相对比例或者滤除某些频率成分的软件或器件。
因此数字滤波器和模拟滤波器有着同样的概念，只是信号形式和实现方式不同。
正因为有该不同点，数字滤波器具有比模拟滤波器精度高、稳定、体积小、重量轻、灵活、不要求阻抗匹配以及实现模拟滤波器无法实现的特殊滤波功能等特点，并且如果要处理的是模拟信号，可通过A/DC和D/AC，在信号形式上进行匹配转换，同样可以使用数字滤波器对模拟信号进行滤波。
数字滤波器实际上是一个线性时不变系统，故可用差分方程、单位脉冲响应h（n）、传输函数（系统函数H（z）及频率响应
来描述。
数字滤波器按照不同分类方法，有许多种类。
但是总体而言，可以分为两类。
一类为经典滤波器，即一般滤波器，其特点是输入信号中有用的频率分量和希望滤波器滤除的频率分量各点有不同的频带，即通过一个合适的选频滤波器达到滤波的目的。
但如果信号与干扰或噪声的频带相互重叠，则显然不能完成有效滤波。
这时，就需要另一类所谓现代滤波器，如维纳滤波器、卡尔曼滤波器、自适应滤波器等最佳地提取信号[1]。
数字滤波器从实现的网络结构或者从单位响应分类，可以分成无限脉冲响应（IIR）滤波器和有限脉冲响应（FIR）滤波器。
它们的系统函数分别为：
（1.1）
（1.2）
（1.1）式中的H（z）称为N阶IIR滤波器的传输函数，（1.2）式为的H（z）称为（N-1）阶FIR滤波器函数。
IIR滤波器产生新的输出，不便需要过去和现在的输入，还需要过去的输出。
而FIR滤波器的输出仅取决于过去的输入，而与过去的输出无关。
IIR滤波器能利用以前所积累的模拟滤波器的成熟的理论及设计图表进行设计的，因而保留了一些典型模拟滤波器优良的幅频特性。
但是设计中只考虑了幅频特性，未考虑相位特性，所设计的滤波器相位特性一般是非线性的。
为了得到线性相位特性，就必须提到FIR滤波器在保证幅度特性满足技术要求的同时，很容易做到有严格的线性相位特性，这也是FIR滤波器的最大的特点。
通常用的数字滤波器一般属于选频滤波器。
假设数字滤波器的传输函数
可用下式来表示：
，其中
称为幅频特性，
称为相频特性。
幅频特性表示信号通过滤波器后各频率分量的衰减情况，而相频特性则反映各频率成分通过滤波器后的相位在时间上的延时情况。
因此，即使两个滤波器幅频特性相同，而相频不一样，对相同的的输入，输出的信号波形也不一样。
一般选频滤波器的技术要求由幅频特性给出，相频特性一般不作要求，但如果对输出波形有要求，则需要考虑相频特性的技术指标，如语音合成、波形传输、图像信号处理等。
如果对输出波形要求有严格要求，则需设计线性相位滤波器。
数字滤波器的设计方法
数字滤波器的设计大致包括三个步骤：
（1）给出所需要的滤波器的技术指标；
（2）设计一
使其逼近所需要的技术指标；
（3）实现所设计的
。
通常，频率选择性滤波器可利用IIR滤波器和FIR滤波器来设计，但是两种不同形式的滤波器的设计不同。
IIR滤波器的单位脉冲响应
为无限长序列，无法由
确定网络结构，而其系统函数
有限，所以设计结果是滤波器系统函数
。
FIR滤波器的设计方法主要是建立在理想滤波器频率特性做某种近似的基础上的，近似方法有窗函数法、频率采样法、优化设计方法等，其设计结果为
。
而对于线型相位滤波器，通常采用FIR滤波器，其单位脉冲响应应满足一定的条件，可以证明其相位特性在整个频带中是严格线性的，这是模拟滤波器无法达到的。
当然，也可以采用IIR滤波器，但必须使用全通网络对其线性相位特性进行相位校正，这样就增加了设计与实现的复杂性。
1.2IIR数字滤波器的设计
IIR数字滤波器的设计步骤如下：
（1）按照一定规则将给出的数字滤波器的技术指标转换成模拟滤波器的技术指标；
（2）根据转换后的技术指标设计模拟低通滤波器
；
（3）再按照一定规则将
转换成
。
如果所设计的数字滤波器是低通型的，那么上述设计工作可以结束，如果所设计的是高通、带通或带阻滤波器，那么还有步骤（4）。
（4）将高通、带通或带阻数字滤波器的技术指标转换成低通模拟滤波器的技术指标，然后按照步骤
（2）设计出低通滤波器
，再将
转换成所需的
。
由上述步骤可知，设计滤波器时，总是先设计模拟低通滤波器，再通过频率变换将之转换成希望设计的滤波器的类型。
模拟滤波器的设计方法已经相当成熟，有着大量的现存图表结果可以查阅，而且MATLAB软件中也包含许多功能强大的设计调用函数，可以用来进行直接设计。
我们可以首先设计一个合适的模拟滤波器，然后变换成满足预定指标的数字滤波器。
这种方法的方便之处在于模拟滤波器已经具有很多简单且现成的设计公式，并且设计参数已经表格化了。
利用模拟滤波器设计数字滤波器，就是要把s平面映射到z平面，这种映射必须满足：
1）s平面的虚轴jΩ必须映射到z平面的单位圆
上；
2）s平面的左半平面Re（s）<0必须映射到z平面的单位圆
的内部︱z︱<1。
一般数字滤波器的设计，我们是先设计几种常用的模拟低通滤波器，高通、带通、带阻等模拟滤波器可利用变量变换方法，由低通滤波器变换得到。
“模拟原型”滤波器有多种设计方法：
如巴特沃斯（Butterworth）滤波器，切比雪夫（Chebyshev）滤波器，椭圆滤波器等。
s-z映射的方法有：
脉冲响应不变法、双线性变换法等。
下表给出了以不同的方法设计IIR滤波器的特性比较。
表1以不同方法设计的IIR滤波器的特性[2]
设计方法
幅度响应
相位响应
滤波器阶数
巴特沃斯
在通带和阻带内的幅度响应平坦，下降慢
近似线性相位
需要较高的阶数
切比雪夫Ⅰ型
过渡速度快于巴特沃斯，但在通带内有更多的纹波
线性介于巴特沃斯与椭圆滤波器之间
阶数介于巴特沃斯与椭圆滤波器之间
切比雪夫Ⅱ型
过渡速度快于巴特沃斯，但在阻带内有更多的纹波
线性介于巴特沃斯与椭圆滤波器之间
阶数介于巴特沃斯与椭圆滤波器之间
椭圆滤波器
在通带和阻带内均为等纹波，
过渡快
高度非线性相位
阶数最低
FIR滤波器的设计问题在于寻求一系统函数
，使其频率响应
逼近滤波器要求的理想频率响应
，其对应的单位脉冲响应
。
1.3用窗函数设计FIR滤波器的基本方法[3]
1.3.1设计思想
从时域从发，设计
逼近理想
。
设理想滤波器
的单位脉冲响应为
。
以低通线性相位FIR数字滤波器为例。
（1.3）
一般是无限长的，且是非因果的，不能直接作为FIR滤波器的单位脉冲响应。
要想得到一个因果的有限长的滤波器h（n），最直接的方法是截断
，即截取为有限长因果序列，并用合适的窗函数进行加权作为FIR滤波器的单位脉冲响应。
按照线性相位滤波器的要求，h（n）必须是偶对称的。
对称中心必须等于滤波器的延时常数，即
（1.4）
用矩形窗设计的FIR低通滤波器，所设计滤波器的幅度函数在通带和阻带都呈现出振荡现象，且最大波纹大约为幅度的9%，这个现象称为吉布斯（Gibbs）效应。
为了消除吉布斯效应，一般采用其他类型的窗函数。
1.3.2．典型的窗函数
（1）矩形窗（RectangleWindow）
其频率响应和幅度响应分别为：
（1.5）
（2）三角形窗（BartlettWindow）
其频率响应为：
（1.6）
（3）汉宁（Hanning）窗，又称升余弦窗
其频率响应和幅度响应分别为：
（1.7）
（4）汉明（Hamming）窗，又称改进的升余弦窗
其幅度响应为：
（1.8）
（5）布莱克曼（Blankman）窗，又称二阶升余弦窗
其幅度响应为：
（1.9）
（6）凯塞（Kaiser）窗
其中：
β是一个可选参数，用来选择主瓣宽度和旁瓣衰减之间的交换关系，一般说来，β越大，过渡带越宽，阻带越小衰减也越大。
I0（·）是第一类修正零阶贝塞尔函数。
若阻带最小衰减表示为
，β的确定可采用下述经验公式：
若滤波器通带和阻带波纹相等即δp=δs时，滤波器节数可通过下式确定：
（1.10）
式中：
1.3.3．设计步骤
利用窗函数设计FIR滤波器的具体步骤如下：
（1）按允许的过渡带宽度△ω及阻带衰减AS，选择合适的窗函数，并估计节数N：
其中A由窗函数的类型决定。
（2）由给定的滤波器的幅频响应参数求出理想的单位脉冲响应
。
（3）确定延时值
（4）计算滤波器的单位取样响应
，
。
（5）验算技术指标是否满足要求。
2快速傅立叶变换（FFT）
2.1FFT算法
FFT，是离散傅立叶变换（DFT）的快速算法，它是根据离散傅立叶变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。
它对傅立叶变换的理论并没有新的发现，但是对于在计算机系统或者说数字系统中应用离散傅立叶变换，可以说是进了一大步。
设有长度为N的序列
的离散傅立叶变换
为：
（2.1）
N点的DFT可以分解为两个N/2点的DFT，每个N/2点的DFT又可以分解为两个N/4点的DFT。
依此类推，当N为2的整数次幂时（
），由于每分解一次降低一阶幂次，所以通过M次的分解，最后全部成为一系列2点DFT运算。
以上就是按时间抽取的快速傅立叶变换（FFT）算法。
当需要进行变换的序列的长度不是2的整数次方的时候，为了使用以2为基的FFT，可以用末尾补零的方法，使其长度延长至2的整数次方。
序列
的离散傅立叶反变换为
（2.2）
离散傅立叶反变换与正变换的区别在于
变为
，并多了一个
的运算。
因为
和
对于推导按时间抽取的快速傅立叶变换算法并无实质性区别，因此可将FFT和快速傅立叶反变换（IFFT）算法合并在同一个程序中。
2．2FFT的优越性
设x（n）为N项的复数序列，由DFT变换，任一X（m）的计算都需要N次复数乘法和N-1次复数加法，而一次复数乘法等于四次实数乘法和两次实数加法，一次复数加法等于两次实数加法，即使把一次复数乘法和一次复数加法定义成一次“运算”（四次实数乘法和四次实数加法），那么求出N项复数序列的X（m）,即N点DFT变换大约就需要N^2次运算。
当N=1024点甚至更多的时候，需要N2=1048576次运算，在FFT中，利用WN的周期性和对称性，把一个N项序列（设N=2k,k为正整数），分为两个N/2项的子序列，每个N/2点DFT变换需要（N/2）^2次运算，再用N次运算把两个N/2点的DFT变换组合成一个N点的DFT变换。
这样变换以后，总的运算次数就变成N+2（N/2）^2=N+N^2/2。
继续上面的例子，N=1024时，总的运算次数就变成了525312次，节省了大约50%的运算量。
而如果我们将这种“一分为二”的思想不断进行下去，直到分成两两一组的DFT运算单元，那么N点的DFT变换就只需要Nlog
（2）（N）次的运算，N在1024点时，运算量仅有10240次，是先前的直接算法的1%，点数越多，运算量的节约就越大，这就是FFT的优越性。
2．3用FFT进行频谱分析
若信号本身是有限长的序列，计算序列的频谱就是直接对序列进行FFT运算求得
，
就代表了序列在
之间的频谱值。
幅度谱
（2.3）
相位谱
（2.4）
若信号是模拟信号，用FFT进行谱分析时，首先必须对信号进行采样，使之变成离散信号，然后就可按照前面的方法用FFT来对连续信号进行谱分析。
按采样定理，采样频率
应大于2倍信号的最高频率，为了满足采样定理，一般在采样之前要设置一个抗混叠低通滤波器。
用FFT对模拟信号进行谱分析的方框图如下所示。
图2.1用FFT进行谱分析方框图
3基于MATLAB的语音信号分析和处理
3.1MATLAB简介
MATLAB是MathWorks公司于1982年推出的一套高性能的数值计算和可视化软件，它集数值分析、矩阵运算、信号处理和图形显示于一体，构成了一个方便的、界面友好的用户环境。
MATLAB的推出得到了各个领域专家学者的广泛关注，其强大的扩展功能为各个领域的应用提供了基础。
借助于MATLAB强大的计算功能进行计算机辅助设计,可以快速有效地设计数字滤波器,大大地简化计算量。
3.2基于MATLAB的语音信号分析
3.2.1语音信号的采集及采样
我们可以利用Windows下的录音机小软件或者其他软件，录制自己一段语音，或者是直接选用一段时间比较短的wav音频文件。
但是这里要注意的是所录制的音频或者是语音要满足一定的要求，如编码格式是PCM编码。
如若编码格式是ADPCM格式（采用魅族MP3自带的外置录制功能录制的语音），则在MATLAB中进行语音读取时会出现错误，提示不能识别。
由于设备有限，故又重新找了一段长为10秒左右的音乐信号进行处理。
其属性如下图所示。
图3.1音乐源文件属性
我们使用MATLAB函数wavread来实现采样。
函数调用格式有多种：
（1）y=wavread（file）,file为所选择的wav文件的路径，返回采样值保存在向量y中；
（2）[y，fs，nbits]=wavread（file），采样值放在向量y中，fs表示采样频率，nbits表示采样位数；
（3）y=wavread（file，N），只读取前N点的采样值放在向量y中；
（4）y=wavread（file,[N1,N2]），读到N1点到N2点的采样值保存在向量y中。
结果如下：
图3.2原始信号波形及频谱图
图3.3用MATLAB得出的fs及bits的值
发现采样频率为音频信号的典型值22050Hz，采样位数与其文件属性的
16位一致，其实在下面的处理中，可以直接将程序中有fs和bits的地方直接用22050及16代替了。
3.2.2语音信号的频谱分析
可以利用MATLAB函数fft对信号进行快速傅里叶变换，得到信号的频谱特性，可以做频谱分析。
建立的M-file如下：
[y,fs,bits]=wavread（'K:
\sound.wav'）;
sound（y,fs,bits）;
Y=fft（y,4096）;
subplot（211）;plot（y）;title（'原始信号的波形'）;
subplot（212）;plot（abs（Y））;title（'原始信号的频谱'）;
以上的sound函数是MATLAB中的一个播放音频的函数，稍后我们可以再重放一下经过滤波后的声音文件，经过比较可以得到经过滤波以后的波形是否有改变。
但是要值得一提是此处的快速傅里叶变换时fft选取了4096点，而采样频率是一秒内采样22050个点。
所以只是选取其中的一部分做频谱分析。
3.2.3设计滤波器进行滤波处理
（1）滤波器的设计
可以用窗函数设计FIR滤波器或者利用双线性变换法等方法设计IIR滤波器。
用matlab设计低通滤波器，首先选用凯塞窗来设计，M-file文件如下：
fp=1000;fc=1200;as=100;ap=1;fs=22050;
wc=2*fc/fs;wp=2*fp/fs;%采用凯塞窗设计
N=ceil（（as-7.95）/（14.36*（wc-wp）/2））+1;
beta=0.1102*（as-8.7）;
win=kaiser（N+1,beta）;
b=fir1（N,wc,win）;
freqz（b,1,512,fs）;
图3.4凯塞窗设计的低通滤波器
图3.5用MATLAB得出凯塞窗设计的滤波器的阶数
从上图可以看出，在低通范围内，滤波器的幅度和线性相位特性满足设计指标，但滤波器的长度，如下图所示，N=708，长度过大，实现起来很困难，主要是滤波器指标太苛刻了。
因此，一般不用窗函数来设计这种类型的滤波器。
其次我们可以用双线性变换法来设计低通滤波器。
用双线性变换法设计的低通滤波器的M-file文件如下：
fp=1000;fc=1200;as=100;ap=1;fs=22050;
wc=2*fc/fs;wp=2*fp/fs;
[n,wn]=ellipord（wp,wc,ap,as）;%此处采用椭圆滤波器原型进行设计
[b,a]=ellip（n,ap,as,wn）;
freqz（b,a,512,fs）;
图3.6用双线性变换法设计的低通滤波器
上面选用椭圆函数设计，滤波器的幅度和相位响应满足设计指标，可以得出滤波器的长度为N=11，实现起来是比较方便的。
（2）用滤波器对采样的音频信号进行滤波
比较以上两种滤波器的性能，用性能比较好的，即N=11的IIR滤波器来对采集的信号进行滤波，滤除高频分量的干扰。
在MATLAB中，FIR滤波器利用函数fftfilt对信号进行滤波，IIR滤波器利用函数filter对信号进行滤波。
具体的实现格式如下：
y=filter（b,a,x）该格式采用数字滤波器对数据进行滤波，既可用于IIR滤波器，也可以用于FIR滤波器。
其中向量b和a分别表示系统函数的分子、分母多项式的系数，若a=1,此时表示FIR滤波器，否则就是IIR滤波器。
该函数是利用给出的向量b和a，对x中的数据进行滤波，结果放入向量y中。
y=fftfilt（b,x）;
该格式是利用基于fft的重叠相加法对数据进行滤波，这种频域滤波技术只对FIR滤波器有效。
该函数能过向量b描述的滤波器对x数据进行滤波。
可以得到滤波后的语音信号的波形和频谱。
（3）滤波程序实现及结果分析
实现滤波的M-file如下：
[y,fs,bits]=wavread（'K:
\sound.wav'）;
Y=fft（y,4096）;
subplot（221）;plot（y）;title（'原始信号的波形'）;
subplot（222）;plot（abs（Y））;title（'原始信号的频谱）;
x=filter（b,a,y）;%滤波的实现
X=fft（x,4096）;
subplot（223）;plot（x）;title（'滤波后信号的波形'）;
subplot（224）;plot（abs（X））;title（'滤波后信号的频谱'）;
sound（x,fs,bits）;%实现滤波后的重放
图3.7滤波前后的信号比较
其中的sound（x,fs,bits）为实现回放语音信号的语句，通过MATLAB处理后可以听到与开始的乐曲相比，声音变得低沉，可见尖锐的高频分量已经被滤除了。
值得说明的是，在MATLAB中使用上述的凯塞窗设计的FIR滤波器来实现滤波时，也可以收到比较好的效果，人耳几乎分辨不出两种滤波的效果，只是在硬件实现时，FIR滤波器由于阶数过高，导致难以实现。
（4）FIR滤波器的滤波分析
以下也给出了使用MATLAB仿真实现的采用上述窗函数实现语音信号滤波及重放的源程序。
[y,fs,bits]=wavread（'K:
\sound.wav'）;
Y=fft（y,4096）;
subplot（221）;plot（y）;title（'原始信号的波形'）;
subplot（222）;plot（abs（Y））;title（'原始信号的频谱）;
fp=1000;fc=1200;as=100;ap=1;fs=22050;
wc=2*fc/fs;wp=2*fp/fs;
N=ceil（（as-7.95）/（14.36*（wc-wp）/2））+1;
beta=0.1102*（as-8.7）;
win=kaiser（N+1,beta）;
b=fir1（N,wc,win）;
x=fftfilt（b,y）;%或者x=filter（b,1,y）;
X=fft（x,4096）;
subplot（223）;plot（x）;title（'滤波后信号的波形'）;
subplot（224）;plot（abs（X））;title（'滤波后信号的频谱'）;
sound（x,fs,bits）;
图3.8FIR滤波器进行滤波前后的信号比较
语音重放后可知效果也不错，而且滤波后的频谱和IIR滤波
升级会员 