音频信号的采样与重构等matlab代码数字信号处理

1、基于MATLAB的音频信号分析与处理 摘要 数字信号处理是一门发展十分迅速、应用非常广泛的前沿学科。它的理论性和实践性都非常强。MATLAB强大的计算仿真功能在数字信号处理领域起着非常大的作用。出于对数字音频处理的兴趣，本文中将尝试利用所学的知识，如采样、滤波、重构等知识，对语音信号或是音频信号进行一定的处理。本文详细给出了利用MATLAB对音频信号进行谱分析，信号滤波和重构的过程，加深了对所学数字信号处理知识的了解。 关键词：滤波 重构 谱分析 Abstract Digital signal processing is an advanced subject which is quickly

2、 developing and widely used. It lays a great emphasis both in theory and practice. MATLAB, the powerful computation and simulation software , plays a great role in digital signal processing field. As for the interest for the digital audio processing, through this paper, I am trying to make some proc

3、essing about the sound (audio) signal with what I have learnt in classroom, such as sampling ,filtering ,reconstruction ,and so on .In this paper ,the processes of frequency amplitude analysis , filtering and construction of the audio signal ,which are based on MATLAB are detailedly presented ,and I

4、 gain more understating about knowledge of digital data processing . Key words : filtering , reconstruction ,frequency amplitude analysis摘要Abstract1 数字滤波器1.1数字滤波器概述1.2 IIR数字滤波器的设计理论1.3 用窗函数设计FIR滤波器1.3.1 设计思想1.3.2典型的窗函数1.3.3设计步骤2 快速傅立叶变换(FFT)2.1 FFT 算法2.2 FFT的优越性 2.3用FFT进行频谱分析3 基于MATLAB的语音信号分析 和处理3.1

5、 MATLAB简介3.2 基于MATLAB的语音信号分析3.2.1 语音信号的采集及采样3.2.2 语音信号的频谱分析3.2.3 设计滤波器进行滤波处理3.3 基于MATLAB的语音信号的谱分析和重构总结参考文献1 数字滤波器1.1 数字滤波器概述数字滤波器是对数字信号实现滤波的线性时不变系统。它通过一定的计算或判断程序来减少干扰信号的比重，它实现的是程序滤波。数字滤波器的输入与输出均为数字信号，通过一定的运算关系改变输入信号所含频率成分的相对比例或者滤除某些频率成分的软件或器件。因此数字滤波器和模拟滤波器有着同样的概念，只是信号形式和实现方式不同。正因为有该不同点，数字滤波器具有 比模拟滤波

6、器精度高、稳定、体积小、重量轻、灵活、不要求阻抗匹配以及实现模拟滤波器无法实现的特殊滤波功能等特点，并且如果要处理的是模拟信号，可通过A/DC和D/AC，在信号形式上进行匹配转换，同样可以使用数字滤波器对模拟信号进行滤波。数字滤波器实际上是一个线性时不变系统，故可用差分方程、单位脉冲响应h(n)、传输函数（系统函数H(z)及频率响应来描述。数字滤波器按照不同分类方法，有许多种类。但是总体而言，可以分为两类。一类为经典滤波器，即一般滤波器，其特点是输入信号中有用的频率分量和希望滤波器滤除的频率分量各点有不同的频带，即通过一个合适的选频滤波器达到滤波的目的。但如果信号与干扰或噪声的频带相互重叠，则

7、显然不能完成有效滤波。这时，就需要另一类所谓现代滤波器，如维纳滤波器、卡尔曼滤波器、自适应滤波器等最佳地提取信号1。数字滤波器从实现的网络结构或者从单位响应分类，可以分成无限脉冲响应(IIR)滤波器和有限脉冲响应(FIR)滤波器。它们的系统函数分别为： (1.1) (1.2) (1.1)式中的H(z)称为N阶IIR滤波器的传输函数，(1.2)式为的H(z)称为（N-1）阶FIR滤波器函数。IIR滤波器产生新的输出，不便需要过去和现在的输入，还需要过去的输出。而FIR滤波器的输出仅取决于过去的输入，而与过去的输出无关。IIR滤波器能利用以前所积累的模拟滤波器的成熟的理论及设计图表进行设计的，因而

8、保留了一些典型模拟滤波器优良的幅频特性。但是设计中只考虑了幅频特性，未考虑相位特性，所设计的滤波器相位特性一般是非线性的。为了得到线性相位特性，就必须提到FIR滤波器在保证幅度特性满足技术要求的同时，很容易做到有严格的线性相位特性，这也是FIR滤波器的最大的特点。通常用的数字滤波器一般属于选频滤波器。假设数字滤波器的传输函数可用下式来表示： ，其中称为幅频特性，称为相频特性。幅频特性表示信号通过滤波器后各频率分量的衰减情况，而相频特性则反映各频率成分通过滤波器后的相位在时间上的延时情况。因此，即使两个滤波器幅频特性相同，而相频不一样，对相同的的输入，输出的信号波形也不一样。一般选频滤波器的技术

9、要求由幅频特性给出，相频特性一般不作要求，但如果对输出波形有要求，则需要考虑相频特性的技术指标，如语音合成、波形传输、图像信号处理等。如果对输出波形要求有严格要求，则需设计线性相位滤波器。数字滤波器的设计方法数字滤波器的设计大致包括三个步骤：(1)给出所需要的滤波器的技术指标；(2)设计一使其逼近所需要的技术指标；(3)实现所设计的。 通常，频率选择性滤波器可利用IIR滤波器和FIR滤波器来设计，但是两种不同形式的滤波器的设计不同。IIR滤波器的单位脉冲响应为无限长序列，无法由确定网络结构，而其系统函数有限，所以设计结果是滤波器系统函数。FIR滤波器的设计方法主要是建立在理想滤波器频率特性做某

10、种近似的基础上的，近似方法有窗函数法、频率采样法、优化设计方法等，其设计结果为 。而对于线型相位滤波器，通常采用FIR滤波器，其单位脉冲响应应满足一定的条件，可以证明其相位特性在整个频带中是严格线性的，这是模拟滤波器无法达到的。当然，也可以采用IIR滤波器，但必须使用全通网络对其线性相位特性进行相位校正，这样就增加了设计与实现的复杂性。1.2 IIR数字滤波器的设计 IIR数字滤波器的设计步骤如下：（1）按照一定规则将给出的数字滤波器的技术指标转换成模拟滤波器的技术指标；（2）根据转换后的技术指标设计模拟低通滤波器；（3）再按照一定规则将转换成。如果所设计的数字滤波器是低通型的，那么上述设计工

11、作可以结束，如果所设计的是高通、带通或带阻滤波器，那么还有步骤（4）。（4）将高通、带通或带阻数字滤波器的技术指标转换成低通模拟滤波器的技术指标，然后按照步骤（2）设计出低通滤波器，再将转换成所需的。由上述步骤可知，设计滤波器时，总是先设计模拟低通滤波器，再通过频率变换将之转换成希望设计的滤波器的类型。模拟滤波器的设计方法已经相当成熟，有着大量的现存图表结果可以查阅，而且MATLAB软件中也包含许多功能强大的设计调用函数，可以用来进行直接设计。我们可以首先设计一个合适的模拟滤波器，然后变换成满足预定指标的数字滤波器。这种方法的方便之处在于模拟滤波器已经具有很多简单且现成的设计公式，并且设计参数

12、已经表格化了。利用模拟滤波器设计数字滤波器，就是要把s平面映射到z平面，这种映射必须满足：1)s平面的虚轴j必须映射到z平面的单位圆上；2)s平面的左半平面Re(s)<0必须映射到z平面的单位圆的内部z <1。一般数字滤波器的设计，我们是先设计几种常用的模拟低通滤波器，高通、带通、带阻等模拟滤波器可利用变量变换方法，由低通滤波器变换得到。“模拟原型”滤波器有多种设计方法：如巴特沃斯(Butterworth)滤波器，切比雪夫(Chebyshev)滤波器，椭圆滤波器等。s-z映射的方法有：脉冲响应不变法、双线性变换法等。下表给出了以不同的方法设计IIR滤波器的特性比较。 表1 以不同方

13、法设计的IIR滤波器的特性2设计方法幅度响应相位响应滤波器阶数巴特沃斯在通带和阻带内的幅度响应平坦，下降慢近似线性相位需要较高的阶数切比雪夫型过渡速度快于巴特沃斯，但在通带内有更多的纹波线性介于巴特沃斯与椭圆滤波器之间阶数介于巴特沃斯与椭圆滤波器之间切比雪夫型过渡速度快于巴特沃斯，但在阻带内有更多的纹波线性介于巴特沃斯与椭圆滤波器之间阶数介于巴特沃斯与椭圆滤波器之间椭圆滤波器在通带和阻带内均为等纹波，过渡快高度非线性相位阶数最低FIR滤波器的设计问题在于寻求一系统函数，使其频率响应逼近滤波器要求的理想频率响应，其对应的单位脉冲响应。1.3 用窗函数设计FIR滤波器的基本方法 3 1.3.1 设

14、计思想从时域从发，设计逼近理想。设理想滤波器的单位脉冲响应为。以低通线性相位FIR数字滤波器为例。 （1.3） 一般是无限长的，且是非因果的，不能直接作为FIR滤波器的单位脉冲响应。要想得到一个因果的有限长的滤波器h(n)，最直接的方法是截断，即截取为有限长因果序列，并用合适的窗函数进行加权作为FIR滤波器的单位脉冲响应。按照线性相位滤波器的要求，h(n)必须是偶对称的。对称中心必须等于滤波器的延时常数，即 (1.4) 用矩形窗设计的FIR低通滤波器，所设计滤波器的幅度函数在通带和阻带都呈现出振荡现象，且最大波纹大约为幅度的9%，这个现象称为吉布斯（Gibbs）效应。为了消除吉布斯效应，一般采

15、用其他类型的窗函数。1.3.2典型的窗函数（1）矩形窗(Rectangle Window) 其频率响应和幅度响应分别为：, (1.5) （2）三角形窗(Bartlett Window) 其频率响应为： (1.6) （3）汉宁(Hanning)窗，又称升余弦窗 其频率响应和幅度响应分别为：(1.7)（4）汉明(Hamming)窗，又称改进的升余弦窗 其幅度响应为： (1.8) （5）布莱克曼(Blankman)窗，又称二阶升余弦窗 其幅度响应为： (1.9) （6）凯塞(Kaiser)窗 其中：是一个可选参数，用来选择主瓣宽度和旁瓣衰减之间的交换关系，一般说来，越大，过渡带越宽，阻带越小衰减也越

16、大。I0(·)是第一类修正零阶贝塞尔函数。 若阻带最小衰减表示为，的确定可采用下述经验公式： 若滤波器通带和阻带波纹相等即p=s时，滤波器节数可通过下式确定： (1.10) 式中：1.3.3设计步骤利用窗函数设计FIR滤波器的具体步骤如下：（1）按允许的过渡带宽度及阻带衰减AS，选择合适的窗函数，并估计节数N：其中A由窗函数的类型决定。（2）由给定的滤波器的幅频响应参数求出理想的单位脉冲响应。（3）确定延时值 （4）计算滤波器的单位取样响应，。（5）验算技术指标是否满足要求。2 快速傅立叶变换(FFT)2.1 FFT 算法FFT，是离散傅立叶变换（DFT）的快速算法，它是根据离散傅立

17、叶变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。它对傅立叶变换的理论并没有新的发现，但是对于在计算机系统或者说数字系统中应用离散傅立叶变换，可以说是进了一大步。 设有长度为N的序列的离散傅立叶变换为： (2.1)N点的DFT可以分解为两个N/2点的DFT，每个N/2点的DFT又可以分解为两个N/4点的DFT。依此类推，当N为2的整数次幂时()，由于每分解一次降低一阶幂次，所以通过M次的分解，最后全部成为一系列2点DFT运算。以上就是按时间抽取的快速傅立叶变换(FFT)算法。当需要进行变换的序列的长度不是2的整数次方的时候，为了使用以2为基的FFT，可以用末尾补零的方法，使其

18、长度延长至2的整数次方。序列的离散傅立叶反变换为 (2.2)离散傅立叶反变换与正变换的区别在于变为，并多了一个的运算。因为和对于推导按时间抽取的快速傅立叶变换算法并无实质性区别，因此可将FFT和快速傅立叶反变换（IFFT）算法合并在同一个程序中。22 FFT的优越性设x(n)为N项的复数序列，由DFT变换，任一X（m）的计算都需要N次复数乘法和N-1次复数加法，而一次复数乘法等于四次实数乘法和两次实数加法，一次复数加法等于两次实数加法，即使把一次复数乘法和一次复数加法定义成一次“运算”（四次实数乘法和四次实数加法），那么求出N项复数序列的X（m）,即N点DFT变换大约就需要N2次运算。当N=1

19、024点甚至更多的时候，需要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+N2/2。继续上面的例子，N=1024时，总的运算次数就变成了525312次，节省了大约50%的运算量。而如果我们将这种“一分为二”的思想不断进行下去，直到分成两两一组的DFT运算单元，那么N点的DFT变换就只需要Nlog(2)(N)次的运算，N在1024点时，运算量仅

20、有10240次，是先前的直接算法的1%，点数越多，运算量的节约就越大，这就是FFT的优越性。23 用FFT进行频谱分析若信号本身是有限长的序列，计算序列的频谱就是直接对序列进行FFT运算求得，就代表了序列在之间的频谱值。幅度谱 （2.3）相位谱 (2.4)若信号是模拟信号，用FFT进行谱分析时，首先必须对信号进行采样，使之变成离散信号，然后就可按照前面的方法用FFT来对连续信号进行谱分析。按采样定理，采样频率应大于2倍信号的最高频率，为了满足采样定理，一般在采样之前要设置一个抗混叠低通滤波器。用FFT对模拟信号进行谱分析的方框图如下所示。抗混叠低通滤波器采样T=1/fsN点FFT图2.1 用F

21、FT进行谱分析方框图 3 基于MATLAB的语音信号分析和处理3.1 MATLAB简介MATLAB是MathWorks公司于1982年推出的一套高性能的数值计算和可视化软件，它集数值分析、矩阵运算、信号处理和图形显示于一体，构成了一个方便的、界面友好的用户环境。MATLAB的推出得到了各个领域专家学者的广泛关注，其强大的扩展功能为各个领域的应用提供了基础。借助于MATLAB 强大的计算功能进行计算机辅助设计,可以快速有效地设计数字滤波器,大大地简化计算量。3.2基于MATLAB的语音信号分析3.2.1 语音信号的采集及采样 我们可以利用Windows 下的录音机小软件或者其他软件，录制自己一段

22、语音，或者是直接选用一段时间比较短的wav音频文件。但是这里要注意的是所录制的音频或者是语音要满足一定的要求 ，如编码格式是PCM编码。如若编码格式是ADPCM格式（采用魅族MP3自带的外置录制功能录制的语音），则在MATLAB中进行语音读取时会出现错误，提示不能识别。由于设备有限，故又重新找了一段长为10秒左右的音乐信号进行处理。其属性如下图所示。 图3.1 音乐源文件属性 我们使用MATLAB函数wavread来实现采样。函数调用格式有多种 ： （1）y=wavread(file),file为所选择的wav文件的路径，返回采样值保存在向量y中；（2） y，fs，nbits=wavread（

23、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对信号进行快速傅里叶变换，得到信号

24、的频谱特性，可以做频谱分析。建立的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点

25、，而采样频率是一秒内采样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);

26、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);

27、%此处采用椭圆滤波器原型进行设计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) 该格式采用数字滤波

28、器对数据进行滤波，既可用于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

29、,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)为实现回放语音

30、信号的语句，通过MATLAB处理后可以听到与开始的乐曲相比，声音变得低沉，可见尖锐的高频分量已经被滤除了。值得说明的是，在MATLAB中使用上述的凯塞窗设计的FIR滤波器来实现滤波时，也可以收到比较好的效果 ，人耳几乎分辨不出两种滤波的效果，只是在硬件实现时，FIR滤波器由于阶数过高，导致难以实现。 （4）FIR滤波器的滤波分析以下也给出了使用MATLAB仿真实现的采用上述窗函数实现语音信号滤波及重放的源程序。 y,fs,bits=wavread('K:sound.wav');Y=fft(y,4096);subplot(221);plot(y);title('原始信号的

31、波形');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('滤波后信号的波形

32、9;);subplot(224);plot(abs(X);title('滤波后信号的频谱');sound(x,fs,bits); 图3.8 FIR滤波器进行滤波前后的信号比较语音重放后可知效果也不错，而且滤波后的频谱和IIR滤波器滤波后的频谱也很类似。3.3 基于MATLAB的语音信号的谱分析和重构还是采用以上的乐曲信号，长度约为10秒，采样率为22050Hz，采样位数为16位。通过数字信号处理的知识可以知道，用fft进行频谱分析时可以得到频谱的特性，从而可以知道信号的能量分布。因为fft频率间隔为2/N，N是离散傅里叶变换用的点数，我们可以利用频谱特性中分析最大值的位置（可以

33、有很多个），它们赌注频率时域的采样有关，相邻的两点之间的距离为f=1/(NTs)，Ts即为采样时间（采样频率的倒数）。前面已经知道了fs=22050Hz。可以能过控制频谱的分量值的大小 ，通过反快速傅里叶变换，重构信号的波形和频谱。其中要用到MATLAB编写程序，要注意是的和C语言中的for语句的不同之处。M-file文件如下： y,fs,bits=wavread('K:sound.WAV');X=0;for i=0:20sample=y(i*4096+1024+1:i*4096+5120);%因为音频文件的长度是10秒，而采样频率为%1024Hz，所以这里采用多段组合的方法，

34、fft点数为4096点%sound(sample,fs,bits);Y=fft(sample,4096);k=1:4096;x(1:4096)=0;%初始化for n=1:4096 if abs(Y(n)>0%此处的语句可以比较精确地控制重构后的效果 P=Y(n); x(n)=P; endendX1=ifft(x,4096); %反快速傅里叶变换进行信号重构X=X,X1; %实现不同分段的信号的拼接endx=fft(X,4096);subplot(311);plot(y(1024+1:20*4096+5120);title('原始信号波形');subplot(312);plot(abs(Y);title('原始信号频谱');subplot(313);plot(X);title('重构后的信号');subplot(414);plot(abs(x);title('重