双音频信号滤波去噪——使用双线性变换法设计的切比雪夫I型滤波器

双音频信号滤波去噪

——使用双线性变换法设计的切比雪夫I型滤波器

摘要本课程设计是采用双线性变换法设计的切比雪夫I型滤波器对双音频信号滤波去噪。在网上下载一段双音频信号，在MATLAB集成环境下，首先用wavread函数求出双音频信号的相关参数，对双音频信号进行读取和加噪；然后再给定相应技术指标，设计一个满足指标的切比雪夫I型滤波器，对该双音频信号进行滤波去噪处理，并绘制对比图，比较滤波前后的波形和频谱并进行分析；最后通过回放双音频信号，对比滤波前后的信号变换。本课程设计成功的对双音频信号进行滤波去噪，初步完成了设计指标。

关键词双音频信号；滤波设计；MATLAB；切比雪夫I型滤波器

1 引言

用麦克风采集一段8000Hz，8k的双音频信号，绘制波形并观察其频谱，给定通带截止频率为2000Hz，阻带截止频率为2150Hz，通带波纹为1dB，阻带波纹为35dB，用双线性变换法设计的一个满足上述指标的切比雪夫I型IIR滤波器，对该双音频信号进行滤波去噪处理。

1.1 课程设计目的

《数字信号处理》课程设计是在学生完成数字信号处理和MATLAB的结合后的基本实验以后开设的。本课程设计的目的是为了让学生综合数字信号处理和MATLAB并实现一个较为完整的小型滤波系统。这一点与验证性的基本实验有本质性的区别。开设课程设计环节的主要目的是通过系统设计、软件仿真、程序安排与调试、写实习报告等步骤，使学生初步掌握工程设计的具体步骤和方法，提高分析问题和解决问题的能力，提高实际应用水平。

1.2 课程设计的要求

（1）学会MATLAB 的使用，掌握MATLAB 的程序设计方法；

（2）滤波器指标必须符合工程实际，根据模拟滤波器的性能指标，确定数字滤波器指标；（3）采用双线性变换法，设计满足上述性能指标要求的ChebyshevI型数字低通滤波器；（4）设计完后应检查其频率响应曲线是否满足指标；

（5）处理结果和分析结论应该一致，而且应符合理论；

（6）独立完成课程设计并按要求编写课程设计报告书；

1.3 设计平台

本次课程设计是在MATLAB软件平台上进行的。MATLAB是矩阵实验室（MATRIX LABORATORY）的简称，是美国MATHWORKS公司推出的具有强大数值分析、矩阵运算、图形绘制和数据处理等功能的软件，现已广泛应用到教学、科研、工程设计等领域[2]。随着MATLAB软件信号处理工具箱的推出，MATLAB已成为信息处理，特别是数字信号处理DSP应用中分析和设计的主要工具。就MATLAB信号处理中的滤波器设计而言，在很大程度上能快速有效地实现滤波器的分析、设计及仿真，大大节约了设计时间，相对传统设计而言，简化了滤波器设计的难度。

2 设计原理

用麦克风采集一段双音频信号，绘制波形并观察其频谱，给定相应技术指标，用双线性变换法设计的切比雪夫I型IIR滤波器，对该双音频信号进行滤波去噪处理，比较滤波前后的波形和频谱并进行分析。

2.1 IIR滤波器

从离散时间来看，若系统的单位抽样(冲激)响应延伸到无穷长，称之为“无限长单位冲激响应系统”，简称为IIR系统。

无限长单位冲激响应（IIR）滤波器有以下几个特点:

(1)系统的单位冲激响应h(n)是无限长；

(2)系统函数H(z)在有限z平面（0

(3)结构上存在着输出到输入的反馈，也就是结构上是递归型的。

IIR滤波器采用递归型结构，即结构上带有反馈环路。同一种系统函数H（z）可以有多种不同的结构，基本网络结构有直接Ⅰ型、直接Ⅱ型、级联型、并联型四种，都具有反馈回路。同时，IIR数字滤波器在设计上可以借助成熟的模拟滤波器的成果，巴特沃斯(Butterworth)滤波器、切比雪夫(Chebyshev)滤波器、椭圆(Cauer)滤波器、贝塞尔(Bessel)滤波器等，这些典型的滤波器各有特点。有现成的设计数据或图表可查，在设计一个IIR数字滤波器时，我们根据指标先写出模拟滤波器的公式，然后通过一定的变换，将模拟滤波器的公式转换成数字滤波器的公式。

2.2 切比雪夫I型滤器

切比雪夫滤波器（又译车比雪夫滤波器）是在通带或阻带上频率响应幅度等波纹波动的滤波器。在通带波动的为“I型切比雪夫滤波器”，在阻带波动的为“II型切比雪夫滤波器”。切比雪夫滤波器在过渡带比巴特沃斯滤波器的衰减快，但频率响应的幅频特性不如后者平坦。切比雪夫滤波器和理想滤波器的频率响应曲线之间的误差最小，但是在通频带内存在幅度波动。

这种滤波器来自切比雪夫多项式，因此得名，用以纪念俄罗斯数学家巴夫尼提·列波维其·切比雪夫

巴特沃兹滤波器在通带内幅度特性是单调下降的，如果阶次一定，则在靠近截止

处，幅度下降很多，或者说，为了使通带内的衰减足够小，需要的阶次N很高，为了克服

这一缺点，采用切比雪夫多项式来逼近所希望的。切比雪夫滤波器的在通带范围内是等幅起伏的，所以在同样的通常内衰减要求下，其阶数较巴特沃兹滤波器要小。切比雪夫滤波器的振幅平方函数为

（2-1）

式中

Ωc—有效通带截止频率

—与通带波纹有关的参量， 大,波纹大 0<

<1

V N （x ）—N 阶切比雪夫多项式

（2-2）

|x|≤1时,|V N (x)|≤1 |x|>1时, |x|↗, V N (x)↗

切比雪夫滤波器的振幅平方特性如图所示,通带内的变化范围为

1(max) →

2

11

ε+

(min)

时，|x|>1,随

↗，

→0 (迅速趋于零)

当

=0时，

（2-3）

N 为偶数，cos 2(

)=1，得到min ，

， （2-4）

N 为奇数，cos 2(

，得到max ，

（2-5）

切比雪夫滤波器的振幅平方特性如图2.1所示。

图2.1 切比雪夫滤波器的振幅平方特性

2.3 双线性变换法

双线性变换法是使数字信号滤波器的频率响应与模拟滤波器的频率响应相似的一种变换方法。为了客服多值映射这一缺点，我们首先把整个s 平面压缩变换到某一中介的s1平面的一条横带里（宽度为

T π2,即从T

π

-到T π),其次再通过上面讨论过的标准变换关系

T s e z 1=将此横带变换到这个z 平面上去，这样就使s 平面与z 平面式一一对应的关系，消

除了多值变换性，也就消除了频谱混叠现象。

将s 平面整个Ωj 轴压缩变换到s1平面Ωj 轴上的T

π-到T π

一段，可以采用以下变换关系：

)2

t a n (1T

Ω=Ω （2-6）

这样，±∞=Ω变换到T

π

±

=Ω1，0=Ω变到01=Ω可将（6）式写成

2

2

221111T j

T j T j

T

j e

e

e e j Ω-ΩΩ-Ω--=

Ω （2-7）

解析延拓到整个s 平面和s1平面，令s j =Ω，11s j =Ω，则得

T

s T

s T j

T j T j

T j e e T s th e

e

e e j 11

111111212

2

22--Ω-ΩΩ-Ω+-=??

????=--=

Ω （2-8） 再将s1平面通过以下标准变换关系映射到z 平面：

T s e z 1= （2-9） 从而得到s 平面和z 平面的单值映射的关系为

1

1

11--+-=z

z s （2-10） s

s

z -+=

11 （2-11） 一般来说，为了使模拟滤波器的某一频率与数字滤波器的任一频率有对应的关系，可以引入待定常数c ，使（6）式和（7）式变换成

)2

t a n (1T

c Ω=Ω （2-12）

T s T s e e c T s c t h s 111121--+-=???

???= （2-13）

仍将T s e z 1=代入（13）式，可得

1

1

11--+-=z

z c s （2-14） s

c s

c z -+=

（2-15） （14）式和（15）式是s 平面与z 平面之间的单值映射关系，这种变换就称为双线性变+换。

3 设计步骤

3.1设计流程图

双音频信号滤波去噪——使用双线性变换法设计的切比雪夫I 型滤波器的设计流程如图3.1所示：

图3.1 双线性变换法切比雪夫I型滤波器对双音频信号去噪流程图

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型滤波器

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); %选择滤波器最小阶数

[b,a] = bilinear(cs,ds,Fs);[C,B,A] = dir2cas(b,a) %双线性变换法实现模拟滤波器到

数字滤波器的转换

[db,mag,pha,grd,w]=freqz_m(b,a); %绘制数字滤波器频率响应幅度图delta=[1,zeros(1,99)];ha=filter(b,a,delta);

得到切比雪夫滤波器的幅度和相位谱如图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函数，对经过切比雪夫I型滤波器之后的信号进行回放，可以听出滤波之后的信号比原始信号更清晰一些，清除了环境噪音。

通过以下语句来进行语音信号回放比较：

sound（x，fs，bits）%播放原始双音频信号

sound（y，fs，bits）%播放加噪后的双音频信号

sound（y_fil，fs，bits）%播放经过滤波处理后的双音频信号

所得结果证明了切比雪夫I型滤波器去噪设计成功。

4 出现的问题及解决方法

在这次的课程设计中，由于理论知识的不踏实以及其他各种原因，我们遇到了不少问题。

（1）在进行双音频信号提取时，进过多次录取才得到理想的双音频信号，在得到理想的波形时，通过多次尝试，和查找书籍及同学讨论，最后猜得到理想的双音频信号的时域图和频谱图

（2）在运用Matlab设计滤波器时，当编辑完前面两条程序时无法放出声音，后来发现我们应当把采集的双音频信号wav文件放到Matlab的work文件夹中。

（3）还要在滤波器性能曲线的As处画一根竖线，这样更方便看出结果。

（4）所有的时间波形横坐标都要化为时间，滤波前后频谱的横坐标应是频率，这样在观察通带截止频率和阻带截止频率时更加精确，误差较小。

5 课程设计心得体会

采用MATLAB设计滤波器，使原来非常繁琐复杂的程序设计变成了简单的函数调用，为滤波器的设和实现开辟了广阔的天地，尤其是MATLAB工具箱使各个领域的研究人员可以直观方便地进行科学研究与工程应用。其中的信号处理工具箱、图像处理工具箱、小波工具箱等更是为数字滤波研究的蓬勃发展提供了可能。MATLAB 信号处理工具箱为滤波器设计及分析提供了非常优秀的辅助设计工具, 在设计数字滤波器时, 善于应用MATLAB进行辅助设计, 能够大大提高设计效率。两周的课程设计让我对MATLAB软件的使用更加的熟练，对切比雪夫滤波器的滤波原理有了更深刻的认识。

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 -1

2.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;