【Matlab】音频信号分析及FIR滤波处理——凯泽(Kaiser)窗

一、前言

1.1 课题内容：

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

1.2 课题要求：

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

1.3 应用价值：

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

二、文献综述

• 数字滤波器有很多种，根据其实现的网络结构或者其冲激响应函数的时域特性，可分为两种：有限冲激响应( FIR，Finite Impulse Response)滤波器和无限冲激响应( IIR，Infinite Impulse Response)滤波器。
• FIR滤波器必须采用间接法，常用的方法有窗函数法、频率采样法和切比雪夫等波纹逼近法。对于线性相位滤波器，经常采用FIR滤波器。在设计FIR滤波器时可以根据对阻带衰减及过渡带的指标要求，选择窗函数类型，并估计窗口长度N，先按照阻带衰减选择窗函数类型。保证阻带衰减满足要求的情况下，尽量选择主瓣的窗函数，再构造希望逼近的频率响应函数，计算。最后加窗便可以得到设计结果。
• FIR滤波器设计方法有窗函数法、频率抽样法等。窗函数设计法比较简单，有闭合形式的公式可循。频率抽样法可以在频域直接设计，在抽样点处与理想滤波器严格相等，可以设计任意幅度响应的滤波器。窗函数设计法是根据设计的性能要求，选择一个理想滤波器，然后用一个合适的窗函数，与理想滤波器在时域中的单位冲激响应相乘，即所谓加窗，得到一个有限长的冲激响应的数字系统，通过调整窗函数的参数来逼近理想滤波器的性能参数，从而达到设计要求。通过窗函数的作用过程可知，一个理想的窗函数在频域的主瓣应该非常窄，有足够的频率分辨率。而旁瓣又能非常低，降低频率之间的干扰。但是在实际中，我们不能同时做到主瓣和旁瓣性能最优，需要在这两者之间取得性能折中。常见的窗函数有矩形(Rectangle)窗、三角形(Fejer)窗、汉宁(Hanning)窗、海明(Hamming)窗、平顶 (Flat Top)窗、凯泽(Kaiser)窗、布莱克曼(Blackman)窗等。矩形窗的优点是主瓣比较集中，频率分辨率最高，缺点是旁瓣较高；三角形窗的主瓣较宽，约等于矩形窗的两倍,但旁瓣小,而且无负旁瓣 ；汉宁窗又称升余弦窗，主瓣变宽，频率分辨率下降，旁瓣减小，有效抑制频谱泄露；海明窗又称为改进升余弦窗，相对汉宁窗来说，其旁瓣更小，但是旁瓣衰减速度变慢；平顶窗在频域通带的波动较小；凯泽窗由一组可调的零阶贝塞尔(Bessel)函数构成，可以通过参数来调整主瓣宽度和旁瓣衰减程度；布莱克曼又称二阶升余弦窗，主瓣宽，旁瓣比较低。
• 窗函数的主瓣宽度和旁瓣峰值衰耗是矛盾的，一项指标的提高总是以另一项指标的下降为代价，窗口选择实际上是对两项指标作权衡。而两项指标是跳变的，于是有人提出可调整窗，适当修改参数，可在这两项指标间作连续的选择。常用的可调整窗是凯泽(Kaiser)窗。凯泽(Kaiser)窗全面地反映主瓣与旁瓣衰减之间的交换关系，可以在它们两者之间自由地选择它们的比重。 表2-1中列出了5种常用的窗函数的特性。
  
• 表2-1中，窗函数在某一个窗长N时，除凯泽(Kaiser)窗以外其他窗函数的系数都是固定的，而凯泽(Kaiser)窗的系数不是固定的，而是随参数值而变化的，凯泽(Kaiser)窗函数的形状也会随着不同的值而变化。
  凯泽(Kaiser)窗在通带波纹和阻带衰减都随参数值而变化，表2-2中列出了部分值与FIR滤波器性能的关系。
  
• 从表2-2中可看出，当参数取不同数值时，阻带的衰减可以从30dB增加到100dB，滤波器的性能与参数的关系极为紧密。所以在滤波器设计中要选择合适的参数值，使滤波器的效果最佳。

三、算法分析

四、算法仿真与结果分析

1. 准备原始音频——录制一段语音信号，时间长度约为6s。将音频信号保存路径为D:\matlab_project\audio.mp3。在MATLAB平台上，用audioread函数调出此语音信号，并得到其音频数据和采样率。

%% 读取音频文件
[x,fs]=audioread('audio.mp3'); % 读取音频信号：x是数据，fs是采样率
x = x(:,1)';                                                  % 取第一列数据转置
N=length(x);                                                  % 计算信号x的长度
t=0:1/fs:(N-1)/fs;                                            % 计算时间范围，样本数除以采样频率
% sound(x,fs);                                                % 播放音频
X=abs(fft(x));X=X(1:N/2);                                     % 对原始信号进行ff变换，截取幅度谱前半部分
deltaf=fs/N;                                                  % 计算频谱的谱线间隔
f=0:deltaf:fs/2-deltaf;                                       % 计算频谱频率范围

在MATLAB中，用plot函数绘制出原始音频信号的时域波形和频谱如下图所示：

2. 对原始音频加入干扰噪声——对原始音频加入窄带高斯噪声。

%% 生成窄带高斯白噪声
noise = wgn(1,N,-20);
[n,Wn,beta,ftype] = kaiserord([10000 10100 10900 11000],...
    [0 1 0],[0.01 0.02 0.01],fs);               % 估计kaiser窗的参数
n=n+rem(n,2);                                   % 保证滤波器系数长N+1为奇数
b_noise=fir1(n,Wn,ftype,kaiser(n+1,beta));      % kaiser窗函数滤波器设计
noise_fil=fftfilt(b_noise,noise);               % 滤出窄带噪声

%% 对原始信号加窄带高斯白噪声
y=x+noise_fil;                   % 将窄带高斯白噪声叠加到原始音频上
% sound(y,fs);                   % 播放加噪信号
Y=abs(fft(y));Y=Y(1:N/2);        % 对加噪信号进行ff变换，截取幅度谱前半部分
figure(2)

在MATLAB中，用plot函数绘制出加噪后的音频信号的时域波形和频谱如图所示

3. 对比信号——对比原始信号与加噪信号频谱图，在10kHz至11kHz有一段窄带信号，则在频率在fc=10000Hz处要滤掉噪声信号，从图中看到，10kHz以后基本没有原始信号的频谱，所以直接设计一个低通滤波器即可。根据频谱图取通带截止频率fp=9900，阻带截止频率fr=10000来设计FIR滤波器。
4. 滤波器设计——要找到一组能满足特定滤波要求的系数向量a和b,而它主要是通过设计指标来实现的。滤波器设计的要求或指标一般是在频域上给出的，常用的滤波器频域指标有：fp、fr、Rp、As。要达到最佳的滤波效果，则需要对fp、fr、Rp、As进行适当的调整。由图可以看出，语音信号可以选择fp=9900;fr=10000;Rp=3;As=50的滤波器。在MATLAB中，通常采用1/2采样频率进行归一化处理。设计程序如程序段4-5所示：

Fs2=fs/2;                                % 奈奎斯特频率
fp=9900;                                 % 通带频率
fr=10000;                                % 阻带频率
Rp=3;                                    % 通带最大波纹衰减
As=40;                                   % 阻带最小衰减
delta1 = (10^(Rp/20)-1)/(10^(Rp/20)+1);  % 通带波纹线性值
delta2 = (1+delta1)*(10^(-As/20));       % 阻带衰减线性值
ff=[fp fr]/Fs2;                          % 频率指标（归一化频率）
A=[1 0];                                 % 和幅值指标
dev=[delta1 delta2];                     % 偏离指标（偏离程度）
[M,Wn,beta,ftype] = kaiserord(ff,A,dev); % 估计kaiser窗的参数
M=M+rem(M,2);                            % 保证滤波器系数长N+1为奇数
n=0:M;                                   % 时间范围
b=fir1(M,Wn,ftype,kaiser(M+1,beta));     % kaiser窗函数滤波器设计
[db,mag,pha,grd,w]=freqz_m(b,1);         % 计算频率响应、相位响应、脉冲响应

运行上面一段程序，通过MATLAB设计数字滤波器基本实现。所设计的滤波器的幅度响应、相位响应和脉冲响应如图所示

由以上图可以看出Rp接近于0，而As大于开始的设定值(As=50db)，由此可见满足了设计要求，达到了设计目的。

5. 音频滤波处理——在MATLAB平台上调用函数fit对信号实现滤波，绘制语音信号去噪前后时域图，频域幅度图形并进行比较。

y_fil=fftfilt(b,y);                        % 用设计好的滤波器对y进行滤波
Y_fil=abs(fft(y_fil)); Y_fil=Y_fil(1:N/2); % 计算频谱取前一半

使用绘图命令得到结果对比图，如图所示

观察上图，滤波后的语音信号发生了衰减，说明滤波器起到了滤波作用，同时通过幅频对比，可以看出滤波器滤掉了一部分频率范围内的信号。最后将滤波后的语音进行回放：播放滤波后的语音信号，发现滤波后的语音信号噪声减小了，同时原始信号强度稍有减弱，没含有尖锐的高频声音，说明滤波已达到设计目的。
6. 结果分析：通过观察滤波前后语音信号波形的变化，原信号的时域图与滤波去噪信号的时域图基本相似，原信号与滤波去噪信号的频谱图波形也大致相似。通过观察可以看到，加噪信号的时域图中大部分都被加入的噪声给遮盖了，加噪信号的频谱图中，我们可以很明显地看到与原信号频谱图相比，它在频率10kHz的地方有一段窄带噪声，而滤波去噪信号的频谱图中该窄带噪声消失，幅度变为了零，波形大致与原图相似，可见滤波去噪效果基本不错。在将三个信号的时域波形和频谱图比较之后，还要通过回放去滤波去噪音频信号，来跟原信号相比，以检验滤波器的效果。在MATLAB中，函数sound函数可以对声音进行回放。回放该滤波去噪信号，便可以感觉到滤波后的语音信号与原信号差不多，该设计基本符合要求。

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

1. 调整通带纹波最大衰减Rp为1dB和阻带最小衰减As为60dB，得到滤波器的响应如图所示：
   

滤波效果如图所示：

可以看到效果和Rp=3和As=50的差别不大。

2. 重新调整通带纹波最大衰减Rp为5dB和阻带最小衰减As为10dB，得到滤波器的响应如图所示：
   
   滤波效果如图所示
   
   效果差的明显，从幅度响应可以看出来通带纹波非常大，最终滤波效果也不好，还有一部分没有被滤除掉。只有在原参数，即Rp=3和As=50的时候滤波效果达到要求，As不能太小。

六、参考文献

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

七、附录

FIR 滤波器——凯泽(Kaiser)窗源程序