【Matlab】音频信号分析及IIR滤波处理——巴特沃斯(Butterworth)

一、前言

1.1 课题内容：

1. 利用麦克风采集语音信号（人的声音、或乐器声乐），人为加上环境噪声(窄带)
2. 分析上述声音信号的频谱，比较两种情况下的差异
3. 根据信号的频谱分布，选取合适的滤波器指标（频率指标、衰减指标），设计对应的 IIR 滤波器
4. 实现数字滤波，将滤波前、后的声音播放，由听觉主观判别滤波效果。并选择、计算合适的客观参数（如：信噪比）评价滤波效果
5. 优化参数，取得更好的滤波效果

1.2 课题要求：

1. 滤波部分要详细讨论各种参数对结果的影响，并对结果数据进行分析、比较与总结
2. 分析、讨论至少两种不同滤波方案的优劣
3. 采用 Matlab 语言编程

1.3 应用价值：

• 20世纪60年代中期形成的一系列数字信号处理方法和算法，如数字滤波器，快速傅里叶变换(FFT)是语音数字信号处理的理论和技术基础。而70年代初期产生的线性预测编码(LPC)算法，为语音信号的数字处理提供了一个强有力的工具。语音信号的编码和压缩是语音信号处理的主要内容。语音信号处理在通信、语音识别与合成、自然语言理解、多媒体数据库以及互联网等多个领域有广泛的应用，同时它对于理解音频类等一般的声音媒体的特点也有很大的帮助。对于移动通信来说，最多的信息是语音信号，语音编码的技术在数字移动通信中具有相当关键的作用，高质量低速率的语音编码技术是数字移动网的永远的追求。所谓语音编码是信源编码，它是将模拟语音信号变成数字信号以便在信道中传输。除了通信带宽的要求外，计算机存储容量的限制也要求对语音信号进行压缩，以满足海量数据情况下进行实时或准实时计算机处理的目的。

二、文献综述

• 数字滤波器有很多种，根据其实现的网络结构或者其冲激响应函数的时域特性，可分为两种：有限冲激响应( FIR，Finite Impulse Response)滤波器和无限冲激响应( IIR，Infinite Impulse Response)滤波器。
• IIR滤波器和FIR滤波器的设计方法不同，IIR滤波器设计方法有间接法和直接法，间接法是借助于模拟滤波器的设计方法进行的。其设计步骤是：先设计过渡模拟滤波器得到系统函数，然后将按某种方法转换成数字滤波器的系统函数。在设计IIR滤波器时可以借助成熟的模拟滤波器的成果，如巴特沃斯、切比雪夫和椭圆滤波器等，根据指标先写出模拟滤波器的公式，然后通过一定的变换，将模拟滤波器的公式转换成数字滤波器的公式。
• IIR滤波器设计方法有巴特沃斯滤波器、切比雪夫滤波器Ⅰ/Ⅱ型、椭圆、贝塞尔滤波器等。巴特沃斯滤波器的特点是通频带内的频率响应曲线最大限度平坦，没有起伏，而在阻频带则逐渐下降为零。在振幅的对数对角频率的波特图上,从某一边界角频率开始,振幅随着角频率的增加而逐步减少,趋向负无穷大。巴特沃斯滤波器的频率特性曲线，无论在通带内还是阻带内都是频率的单调函数。因此，当通带的边界处满足指标要求时，通带内肯定会有裕量。所以，更有效的设计方法应该是将精确度均匀的分布在整个通带或阻带内，或者同时分布在两者之内。这样就可用较低阶数的系统满足要求。巴特沃斯滤波器通带内的幅频响应曲线能得到最大限度的平滑，但牺牲了截止频率的坡度。切比雪夫滤波器是在通带或阻带上频率响应幅度等波纹波动的滤波器，振幅特性在通带内是等波纹。在阻带内是单调的称为切比雪夫I型滤波器；振幅特性在通带内是单调的，在阻带内是等波纹的称为切比雪夫II型滤波器。采用何种形式的切比雪夫滤波器取决于实际用途。切比雪夫I型滤波器通带内等波纹，阻带内单调；切比雪夫Ⅱ型滤波器通带内单调，然而阻带内等波纹。椭圆滤波器阻带和通带内都是等波纹的，但下降的坡度更大，而且可以以更低的阶数实现和其他两类滤波器一样的性能指标。

三、算法分析

• IIR滤波器与FIR滤波器相比，具有相位特性差的缺点，但它的的结构简单、运算量小，具有经济、高效的特点，并且可以用较少的阶数获得很高的选择性。因此也得到了广泛应用。IIR滤波器具有无限长度的单位脉冲响应，在结果上存在反馈回路，具有递归性，即IIR滤波器的输出不仅与输入有关，而且与过去的输出有关。巴特沃斯滤波器相较于切比雪夫滤波器、椭圆、贝塞尔等经典滤波器，在同等阶数条件下有通频带内的频率响应曲线最大限度平坦，没有起伏，而在阻频带则逐渐下降为零的特点。
• 目前IIR数字滤波器设计的最通用方法是借助于模拟滤波器的设计方法。模拟滤波器已经有了一套相当成熟的方法，它不但有完整的设计公式，而且还有较为完整的图表供查询，因此，充分利用这些已经有的资源会给数字滤波器的设计代理很大的方便。
  按照这一方式,巴特沃斯IIR数字滤波器的设计过程为:

1. 确定数字滤波器实际需满足的技术指标参数，包括通带截至频率Ω_p、阻带起始频率Ω_s、通带最大波动A_p、阻带最小衰减A_s
2. 计算阶数N
   
3. 根据阶数查询表，得到归一化系统函数H_α(P)的分母系数
   
   
4. 去归一化得到H_α(s)，至此得到模拟巴特沃斯滤波器。
5. 利用双线性变换法，将s平面映射到z平面，得到数字巴特沃斯IIR滤波器。

四、算法仿真与结果分析

1. 使用MATLAB软件中recordblocking函数录制一段4秒钟，采样频率为48000HZ的音频，录制完成后存为“Voice.wav”放置于桌面上。再使用同样的参数，在录制音频的同时加入背景噪声，录制完成后保存为“VoiceAddNoise.wav”放置于桌面上。
2. 在MATLAB中使用以下程序读取音频信号

1.[voice_noise,Fs] = audioread('VoiceAddNoise.wav');  %音频文件
2.[voice,Fs] = audioread('Voice.wav');  %音频文件
3.t=0:1/Fs:1;                 % 时间轴 
4.len=length(voice_noise);
5.y1=fft(voice,len);          % 对信号做len点FFT
6.y2=fft(voice_noise,len);   % 对信号做len点FFT
7.f=Fs*(0:len/2-1)/len;

3. 对其进行快速傅里叶变换，如下所示
   
4. 对比“Voice”与“VoiceAddNoise”的幅度谱图，发现语音信号频率集中在0-3kHZ的范围内，则在频率大于3000Hz的是需要滤掉噪声信号，于是设置通带截止频率为3000HZ，阻带截止频率为3800HZ，通带纹波最大衰减为3dB，阻带衰减为30dB通过以上参数设计IIR滤波器。
5. 计算最小滤波器阶数N，调用MATLAB中butter函数快速设计低通滤波器，如下所示。

1.% 3.IIR滤波器设计
2.N=0;                                  % 阶数
3.Fp=3000;                              % 通带截止频率1000Hz
4.Fc=3800;                              % 阻带截止频率2000Hz
5.Rp=3;                                 % 通带波纹最大衰减为1dB
6.Rs=30;                                % 阻带衰减为60dB
7.% 3.0计算最小滤波器阶数
8.na=sqrt(10^(0.1*Rp)-1);             % 开根号
9.ea=sqrt(10^(0.1*Rs)-1);             % 开根号
10.N=ceil(log10(ea/na)/log10(Fc/Fp)); % 朝正无穷大方向取整
11.% 3.1巴特沃斯滤波器
12.Wn=Fp*2/Fs;
13.[Bb,Ba]=butter(N,Wn,'low');		 % 调用MATLAB_butter函数快速设计滤波器
14.[BH,BW]=freqz(Bb,Ba);               % 绘制频率响应曲线
15.Bf=filter(Bb,Ba,voice_noise);      % 进行低通滤波

6. 滤波器设计完成后，在MATLAB平台上调用函数filter对信号实现低通滤波，绘制语音信号去噪前后时域图，频域幅度图形并进行比较。通过频谱图可以看出，高于3kHZ的信号被削减，同时播放去噪后的语音文件，发现蝉鸣声消失，只留下人声信号。
   
   

五、关键参数的调试与分析

1. 设置通带截至频率为3300HZ，阻带截至频率为3800HZ，通带纹波最大衰减3dB,阻带衰减30dB,由频谱图看出在0-3kHz的语音信号频谱中出现了失真现象，如下图所示。
   
2. 调整通带截止频率为3000HZ，阻带截止频率为3800HZ，通带纹波最大衰减为3dB，阻带衰减为30dB，尽管有部分语音信号被削弱，但整体未出现大幅度失真，并且成功滤除高于3kHZ的蝉鸣信号，如下图所示。
   

六、参考文献

[1] 庞建丽，高丽娜.基于Matlab的IIR数字滤波器设计方法比较及应用[J].黄淮学院(现代电子技术)，2010 (11): 103-105
[2] 岳学东，买鹏.基于MATLAB的FIR数字滤波器最优设计[J].中国科技信息，2005
[3] 程佩青.数字信号处理教程[M].北京：清华大学出版社，2002
[4] 赵晓群，张洁.巴特沃斯低通滤波器的实现方法研究[J]。大连民族学院学报，2013,15(01):72-75.

七、附录

IIR滤波器——巴特沃斯(Butterworth) 源程序