语音信号滤波去噪
——窗函数法设计FIR低通滤波器
摘 要 本课程设计是采用窗函数法设计FIR低通滤波器对双音频信号滤波去噪。录制一段语音信号，在MATLAB集成环境下，首先用wavread函数求出语音信号的相关参数，对语音信号进行读取和加噪；然后再给定相应技术指标，设计一个满足指标的FIR低通滤波器，对该双音频信号进行滤波去噪处理，并绘制对比图，比较滤波前后的波形和频谱并进行分析；最后通过回放双音频信号，对比滤波前后的信号变换。本课程设计成功的对双音频信号进行滤波去噪，初步完成了设计指标。
关键词 双音频信号； 滤波设计； MATLAB；FIR低通滤波器
1 引 言
用麦克风采集一段8000Hz，8k的双音频信号，绘制波形并观察其频谱，给定通带截止频率为2000Hz，阻带截止频率为2150Hz，通带波纹为1dB，阻带波纹为35dB，用双线性变换法设计的一个满足上述指标的切比雪夫I型IIR滤波器，对该双音频信号进行滤波去噪处理。
1.1 课程设计目的
《数字信号处理》课程设计是在学生完成数字信号处理和MATLAB的结合后的基本实验以后开设的。本课程设计的目的是为了让我们综合数字信号处理和MATLAB并实现一个较为完整的小型滤波系统。这一点与验证性的基本实验有本质性的区别。开设课程设计环节的主要目的是通过系统设计、软件仿真、程序安排与调试、写实习报告等步骤，使学生初步掌握工程设计的具体步骤和方法，提高分析问题和解决问题的能力，提高实际应用水平。
1.2 课程设计的要求
（1）学会 MATLAB 的使用，掌握 MATLAB 的程序设计方法；
（2）滤波器指标必须符合工程实际，根据模拟滤波器的性能指标，确定数字滤波器指标；
（3）采用窗函数法，设计满足上述性能指标要求的FIR型数字低通滤波器；
（4）设计完后应检查其频率响应曲线是否满足指标；
（5）处理结果和分析结论应该一致，而且应符合理论；
（6）独立完成课程设计并按要求编写课程设计报告书；
1.3 设计平台
本次课程设计是在MATLAB软件平台上进行的。MATLAB是矩阵实验室（MATRIX LABORATORY）的简称，是美国MATHWORKS公司推出的具有强大数值分析、矩阵运算、图形绘制和数据处理等功能的软件，现已广泛应用到教学、科研、工程设计等领域[2]。随着MATLAB软件信号处理工具箱的推出，MATLAB已成为信息处理，特别是数字信号处理DSP应用中分析和设计的主要工具。就MATLAB信号处理中的滤波器设计而言，在很大程度上能快速有效地实现滤波器的分析、设计及仿真，大大节约了设计时间，相对传统设计而言，简化了滤波器设计的难度。
2 设计原理
2.1用麦克风采集一段双音频信号，绘制波形并观察其频谱，给定相应技术指标，用双线性变换法设计的窗函数法设计FIR低通滤波器，对该双音频信号进行滤波去噪处理，比较滤波前后的波形和频谱并进行分析。
2.2 模拟含噪语音信号合成；
在MATLAB软件平台下，给原始的语音信号叠加上噪声，噪声类型分为如下几种：a白噪声；b单频噪色（正弦干扰）；c多频噪声（多正弦干扰）；d其它干扰。绘出叠加噪声后的语音信号时域和频谱图，在视觉上与原始语音信号图形对比，也可通过Windows播放软件从听觉上进行对比，分析含噪语音信号频谱和时域波形的改变。
2.3 数字滤波器设计及滤波；
已知数字滤波器的性能指标为：，采用窗函数法设计低通型FIR滤波器来对叠加噪声前后的语音信号进行滤波处理，绘出滤波器的频域响应、滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比，分析信号的变化。分别绘制出理想冲激响应和实际冲激响应结果图。并且给出幅度响应结果图。
3 设计步骤
3.1设计流程图
语音信号滤波去噪——使用窗函数法设计FIR低通滤波器的设计流程如图3.1所示：
图3.1 窗函数法设计FIR低通滤波器对双音频信号去噪流程图
3.2语音信号的采集
点击windows系统桌面的“开始”按钮，点击开始菜单栏里的“附件”，选择“录音机”选项，点击录音机“文件”选项，进入“声音选定”设置，把属性一栏设置成“8000Hz，8位，单声道，7KB/秒”(见图3.2)。点击确定，然后开始双音频信号的采集，采集时间为2—3秒左右为最佳。采集的声音文件以“.wav”格式存储（见图3.3）。
图3.2 采集声音的参数设置
图3.3 采集声音
3.3语音信号的频谱分析
在MATLAB中编辑m函数，使用wavread函数读取采集的声音文件（.wav)将它赋值给某一向量，再对其进行采样,然后使用plot语句画出相关的频谱图形。
（1）Wavread函数调用格式：
[y,fs,nbits]=wavread(file)
功能说明：采样值放在向量x中，fs表示采样频率（Hz），bits表示采样位数。
（2）快速傅里叶变换算法FFT计算DFT的函数fft,其调用格式如下：
Xk=fft(x,n)
参数x为被变换的时域序列向量，N是DFT变换区间长度，当n大于x的长度时，fft函数自动在x后面补零。，当n小于xn的长度时，fft函数计算x的前n个元素，忽略其后面的元素。
在本次课程设计中，我们利用fft函数对双音频信号进行快速傅里叶变换，就可以得到信号的频谱特性。
（3）声音采样文件读取的程序（文件名：xiaomiao.wav)；
双音频信号的提取：
[x,fs,bits]=wavread(‘D:\MATLAB7\work\xiaomiao.wav’); % 输入参数为文件的全路径和文件名，输出的第一个参数是每个样本的值，fs是生成该波形文件时的采样率，bits是波形文件每样本的编码位数。
sound(x,fs,bits); % 按指定的采样率和每样本编码位数回放
N=length(x); % 计算信号x的长度
fn=2200; % 单频噪声频率，此参数可改
t=0:1/fs:(N-1)/fs; % 计算时间范围，样本数除以采样频率
x=x(:,1)’; y=x+sin(fn*2*pi*t);
sound(y,fs,bits); % 应该可以明显听出有尖锐的单频啸叫声
X=abs(fft(x)); Y=abs(fft(y)); % 对原始信号和加噪信号进行fft变换,取幅度谱
X=X(1:N/2); Y=Y(1:N/2); % 截取前半部分
deltaf=fs/N; % 计算频谱的谱线间隔
f=0:deltaf:fs/2-deltaf; % 计算频谱频率范围
得到的原始语音信号和加上噪音后的语音信号的时域波形和频谱图如图3.4所示。
图3.4 原始双音频信号和加噪后信号的时域和频谱图
3.4滤波器设计
设计指标：通带截止频率为2000Hz，阻带截止频率为2150Hz，通带波纹为1dB，阻带波纹为35dB，用双线性变换法设计的一个满足上述指标的切比雪夫I型滤波器
双线性变换法设计切比雪夫I型滤波器
fs=22050;
[x,fs,bits]=wavread(‘D:\我的祖国22.wav’);
sound(x,fs,bits); % 按指定的采样率和每样本编码位数回放
fn=2200; % 单频噪声频率，此参数可改
t=0:1/fs:(N-1)/fs; % 计算时间范围，样本数除以采样频率
%x=x(:,1)’;
y=x+awgn(x,10);
N=length(y);
%sound(y,fs,bits); % 应该可以明显听出有尖锐的单频啸叫声
X=abs(fft(x)); Y=abs(fft(y)); % 对原始信号和加噪信号进行fft变换,取幅度谱
%X=X(1:N/2); Y=Y(1:N/2); % 截取前半部分
%deltaf=fs/N; % 计算频谱的谱线间隔
%f=0:deltaf:fs/2-deltaf; % 计算频谱频率范围
y_p=fft(y,N);
f=fs*(0:N/2-1)/N;
figure(1);
subplot(2,1,1);
plot(y);
title(‘原始语音信号采样后时域波形’);
xlabel(‘时间轴’);
ylabel(‘幅度A’);
subplot(2,1,2);
plot(f,abs(y_p(1:N/2)));
title(‘原始语音信号采样后频谱图’);
xlabel(‘频率HZ’);
ylabel(‘频率幅值’);
得到窗函数法设计FIR低通滤波器的幅度和相位谱如图3.5所示。
图3.5 切比雪夫滤波器的幅度和相位谱
3.5双音频的滤波程序及滤波效果图
滤波程序：
y_fil=filter(b,a,y); % 用设计好的滤波器对y进行滤波
Y_fil=abs(fft(y_fil));Y_fil=Y_fil(1:N/2); % 计算频谱取前一半
画出滤波前后双音频信号的时域和频谱图如图3.6所示。
图3.6 滤波前后双音频信号的时域和频谱图
3.6结果分析
在MATLAB中，经sound函数，对经过窗函数法设计FIR低通滤波器之后的信号进行回放，可以听出滤波之后的信号比原始信号更清晰一些，清除了环境噪音。
通过以下语句来进行语音信号回放比较：
sound（x，fs，bits） %播放原始双音频信号
sound（y，fs，bits） %播放加噪后的双音频信号
sound（y_fil，fs，bits） %播放经过滤波处理后的双音频信号
所得结果证明了窗函数法设计FIR低通滤波器去噪设计成功。
4 出现的问题及解决方法
在这次的课程设计中，由于理论知识的不踏实以及其他各种原因，我们遇到了不少问题。
（1）在进行双音频信号提取时，进过多次录取才得到理想的双音频信号，在得到理想的波形时，通过多次尝试，和查找书籍及同学讨论，最后猜得到理想的双音频信号的时域图和频谱图
（2）在运用Matlab设计滤波器时，当编辑完前面两条程序时无法放出声音，后来发现我们应当把采集的双音频信号wav文件放到Matlab的work文件夹中。
（3）还要在滤波器性能曲线的As处画一根竖线，这样更方便看出结果。
（4）所有的时间波形横坐标都要化为时间，滤波前后频谱的横坐标应是频率，这样在观察通带截止频率和阻带截止频率时更加精确，误差较小。
5 课程设计心得体会
采用MATLAB设计滤波器，使原来非常繁琐复杂的程序设计变成了简单的函数调用，为滤波器的设和实现开辟了广阔的天地，尤其是MATLAB工具箱使各个领域的研究人员可以直观方便地进行科学研究与工程应用。其中的信号处理工具箱、图像处理工具箱、小波工具箱等更是为数字滤波研究的蓬勃发展提供了可能。MATLAB 信号处理工具箱为滤波器设计及分析提供了非常优秀的辅助设计工具, 在设计数字滤波器时, 善于应用MATLAB进行辅助设计, 能够大大提高设计效率。两周的课程设计让我对MATLAB软件的使用更加的熟练，对窗函数法设计FIR低通滤波器原理有了更深刻的认识。
6 参考文献
[1] 程佩青．数字信号处理教程［M］．北京：清华大学出版社，2002
[2] 薛年喜主编．MATLAB 在数字信号处理中的应用[M]．北京：清华大学出版社，2002
[3] 维纳•K•恩格尔，约翰•G•普罗克斯.《数字信号处理》[M].西安交通大学出版社，2002
[4] 董长虹等. MATLAB信号处理与应用[M].北京：国防工业出版社，2005
[5] [美] M.H.海因斯 著，张建华等 译.数字信号处理[M].北京：科学出版社，2002
[6] 张葛祥，李 娜. MATLAB仿真技术与应用[M].北京：清华大学出版社，2007
附录一: 源程序设计清单
原始双音频信号的采集及加噪：
[x,fs,bits]=wavread(‘D:\MATLAB7\work\xiaomiao.wav’); % 输入参数为文件的全路径和文件名，输出的第一个参数是每个样本的值，fs是生成该波形文件时的采样率，bits是波形文件每样本的编码位数。
sound(x,fs,bits); % 按指定的采样率和每样本编码位数回放
N=length(x); % 计算信号x的长度
fn=2200; % 单频噪声频率，此参数可改
t=0:1/fs:(N-1)/fs; % 计算时间范围，样本数除以采样频率
x=x(:,1)’; y=x+sin(fn*2*pi*t);
sound(y,fs,bits); % 应该可以明显听出有尖锐的单频啸叫声
X=abs(fft(x)); Y=abs(fft(y)); % 对原始信号和加噪信号进行fft变换,取幅度谱
X=X(1:N/2); Y=Y(1:N/2); % 截取前半部分
Warning: Integer operands are required for colon operator when used as index.
Warning: Integer operands are required for colon operator when used as index.
deltaf=fs/N; % 计算频谱的谱线间隔
f=0:deltaf:fs/2-deltaf; % 计算频谱频率范围
subplot(2,2,1);plot(t,x);xlabel(‘时间(单位:s)’);ylabel(‘幅度’);title(‘原始双音频信号’);grid
subplot(2,2,2);plot(f,X);xlabel(‘频率(单位:Hz)’);ylabel(‘幅度谱’);title(‘原始双音频信号幅度谱图’); axis([0,6*10^3,0,6000]);grid
subplot(2,2,3);plot(t,y);xlabel(‘时间(单位:s)’);ylabel(‘幅度’);title(‘加噪后的双音频信号’);grid
subplot(2,2,4);plot(f,Y);xlabel(‘频率(单位:Hz)’);ylabel(‘幅度谱’);title(‘加噪后的双音频信号幅度谱图’); axis([0,6*10^3,0,6000]);grid
使用双线性变换法设计切比雪夫I型滤波器：
fp=fn-200;fc=fn-50; %定义通带和阻带截止频率
Rp=1;As=35; % 定义通带波纹和阻带衰减
wp=fp/fs*2*pi; ws=fc/fs*2*pi; %计算对应的数字频率
T=1;Fs = 1/T; %定义采样间隔
Omegap=(2/T)*tan(wp/2);
Omegas=(2/T)*tan(ws/2); %截止频率预畸变
[cs,ds] = afd_chb1(Omegap,Omegas,Rp,As); %选择滤波器最小阶数
Warning: Function call afd_chb1 invokes inexact match d:\MATLAB7\work\AFD_CHB1.M.
*** 切比雪夫-1 滤波器阶次 = 14
Warning: Function call u_chb1ap invokes inexact match d:\MATLAB7\work\U_CHB1AP.M.
> In AFD_CHB1 at 29
[b,a] = bilinear(cs,ds,Fs);[C,B,A] = dir2cas(b,a) %双线性变换法实现模拟滤波器到
数字滤波器的转换
Warning: Function call dir2cas invokes inexact match d:\MATLAB7\work\DIR2CAS.M.
C =
3.3307e-016
B =
1.0000 -0.0321 0.4021
1.0000 -0.2091 0.0238
1.0000 -0.2206 5.8342
1.0000 -0.3400 0.1487
1.0000 -0.5438 0.7109
1.0000 -2.1153 8.5100
1.0000 -12.5391 364.1406
A =
1.0000 -1.9050 0.9776
1.0000 -1.9116 0.9540
1.0000 -1.9147 0.9950
1.0000 -1.9185 0.9754
1.0000 -1.9315 0.9425
1.0000 -1.9354 0.9570
1.0000 -1.9430 0.9444
[db,mag,pha,grd,w]=freqz_m(b,a); %绘制数字滤波器频率响应幅度图
Warning: Function call freqz_m invokes inexact match d:\MATLAB7\work\FREQZ_M.M. delta=[1,zeros(1,99)];ha=filter(b,a,delta);
figure
Subplot(221);plot(w/pi,db); title(‘切比雪夫I 滤波器幅度(dB)’);
xlabel(‘w’);ylabel(‘dB’);axis([0,1,-150,10]);grid
Subplot(222);plot(w/pi,mag);title(‘切比雪夫I滤波器幅度响应’);xlabel(‘w’);ylabel(‘幅值 |H|’); axis([0,1,0,1]);grid
Subplot(223);plot(w/pi,pha);title(‘切比雪夫I 滤波器相位响应’); xlabel(‘w’);ylabel(‘相位（单位：π ）’);axis([0,1,-4,4]);grid
对加噪信号进行滤波处理并画出加噪前后的信号时域和频谱：
y_fil=filter(b,a,y); % 用设计好的滤波器对y进行滤波
Y_fil=abs(fft(y_fil));Y_fil=Y_fil(1:N/2); % 计算频谱取前一半
Warning: Integer operands are required for colon operator when used as index.
figure
subplot(3,2,1);plot(t,x);xlabel(‘时间(单位:s)’);ylabel(‘幅度’);title(‘原始双音频信号’);grid
subplot(3,2,2);plot(f,X);xlabel(‘频率(单位:Hz)’);ylabel(‘幅度谱’);title(‘原始双音频信号幅度谱图’);grid
axis([0,6*10^3,0,6000])
subplot(3,2,3);plot(t,y);xlabel(‘时间(单位:s)’);ylabel(‘幅度’);title(‘加噪后的双音频信号’);grid
subplot(3,2,4);plot(f,Y);xlabel(‘频率(单位:Hz)’);ylabel(‘幅度谱’);title(‘加噪后的双音频信号幅度谱图’);grid
axis([0,6*10^3,0,6000])
subplot(3,2,5);plot(t,y_fil);xlabel(‘时间(单位:s)’);ylabel(‘幅度’);title(‘滤波后的双音频信号’);grid
axis([0,6,-1,1])
subplot(3,2,6);plot(f,Y_fil);xlabel(‘频率(单位:Hz)’);ylabel(‘幅度谱’);title(‘滤波后的双音频信号幅度谱图’);grid
axis([0,6*10^3,0,6000])
附录二：afd_chb1函数
function [b,a] = afd_chb1(Wp,Ws,Rp,As);
% 切比雪夫-1型模拟低通滤波器设计
% —————————————–
% [b,a] = afd_chb1(Wp,Ws,Rp,As);
% b = Ha(s) 分子的系数
% a = Ha(s) 分母的系数
% Wp = 以弧度/秒为单位的通带边缘频率; Wp > 0
% Ws = 以弧度/秒为单位的阻带边缘频率; Ws > Wp > 0
% Rp = 通带中的振幅波动的+dB数; (Rp > 0)
% As = 阻带衰减的+dB数; (As > 0)
%
if Wp <= 0
error(‘通带边缘必须大于 0’)
end
if Ws <= Wp
error(‘阻带边缘必须大于通带边缘’)
end
if (Rp <= 0) | (As < 0)
error(‘通带波动或阻带衰减必须大于0’)
end
ep = sqrt(10^(Rp/10)-1);
A = 10^(As/20);
OmegaC = Wp;
OmegaR = Ws/Wp;
g = sqrt(A*A-1)/ep;
N = ceil(log10(g+sqrt(g*g-1))/log10(OmegaR+sqrt(OmegaR*OmegaR-1)));
fprintf(‘\n*** 切比雪夫-1 滤波器阶次 = %2.0f \n’,N)
[b,a]=u_chb1ap(N,Rp,OmegaC);
附录三：freqz.m.m函数
function [db,mag,pha,grd,w] = freqz_m(b,a);
% freqz 子程序的改进版本
% ————————————
% [db,mag,pha,grd,w] = freqz_m(b,a);
% db = [0 到pi弧度]区间内的相对振幅(db)
% mag = [0 到pi弧度]区间内的绝对振幅
% pha = [0 到pi弧度]区间内的相位响应
% grd = [0 到pi弧度]区间内的群迟延
% w = [0 到pi弧度]区间内的501个频率样本向量
% b = Ha(z)的分子多项式系数(对FIR b=h)
% a = Ha(z)的分母多项式系数(对 FIR: a=[1])
%
[H,w] = freqz(b,a,1000,’whole’);
H = (H(1:1:501))’; w = (w(1:1:501))’;
mag = abs(H);
db = 20*log10((mag+eps)/max(mag));
pha = angle(H);
% pha = unwrap(angle(H));
grd = grpdelay(b,a,w);
% grd = diff(pha);
% grd = [grd(1) grd];
% grd = [0 grd(1:1:500); grd; grd(2:1:501) 0];
% grd = median(grd)*500/pi;