基于MATLAB的数字信号滤波器设计

一、课题目的

1、学会MATLAB的使用,掌握其程序设计方法,学会对信号进行分析和处理;

2、掌握语音信号的采集、存储和时频分析;

3、要求掌握IIR数字滤波器的设计原理、设计方法和设计步骤;

4、学习用窗函数法设计FIR数字滤波器的原理及其设计步骤;

5、了解学习GUI图形用户界面,进行实验程序的演示。

二.实验内容

1、录制一段个人自己的语音信号,并对录制的信号进行采样;画出采样后语音信号的时域波形和频谱图;

2、给定滤波器的性能指标,采用窗函数法和双线性变换法设计滤波器,并画出滤波器的频率响应;

3、用设计的滤波器对采集的信号进行滤波,画出滤波后信号的时域波形和频谱,并对滤波前后的信号进行对比,分析信号的变化;回放语音信号;

4、换一个性别相异的人录制同样一段语音内容,分析两段内容相同的语音信号频谱之间的特点;

5、录制一段同样长时间的背景噪声叠加到语音信号中,分析叠加前后信号频谱的变化,设计一个合适的滤波器,能够把该噪声滤除。

三.实验原理

分析实验的几个要求可以看出,实验主要考察的是数字滤波器的设计及语音信号的滤波,涉及到男声和女声的分析。从本质上说,本实验涉及的信号是人的语言,众所周知,人声的频率范围是20Hz到20000Hz,进一步的分析知声音范围大约在65-950 Hz,所以设计的是低通滤波器。

1.语音信号的采集

熟悉并掌握MATLAB中有关声音(wave)录制、播放、存储和读取的函数,在MATLAB环境中,有关声音的函数有:

a:y=wavrecord(N,fs,Dtype);利用系统音频输入设备录音,以fs为采样频率,默认值为11025,即以11025HZ进行采样。Dtype为采样数据的存储格式,用字符串指定,可以是:‘double’、‘single’、‘int16’、‘int8’其中只有int8是采用8位精度进行采样,其它三种都是16位采样结果转换为指定的MATLAB数据;

b:wavplay(y,fs);利用系统音频输出设备播放,以fs为播放频率,播放语音信号y;

c:wavwrite((y,fs,wavfile);创建音频文件;

d:y=wavread(file);读取音频文件;

关于声音的函数还有sound();soundsc();等。

2.滤波器的设计原理:

数字滤波器是数字信号处理中及其重要的一部分。随着信息时代和数字技术的发展,受到人们越来越多的重视。数字滤波器可以通过数值运算实现滤波,所以数字滤波器处理精度高、稳定、体积小、重量轻、灵活不存在阻抗匹配问题,可以实现模拟滤波器无法实现的特殊功能。数字滤波器种类很多,根据其实现的网络结构或者其冲激响应函数的时域特性,可分为两种,即有限冲激响应( FIR,Finite Impulse Response)滤波器和无限冲激响应( IIR,Infinite Impulse Response)滤波器。

FIR滤波器结构上主要是非递归结构,没有输出到输入的反馈,系统函数H (z)在处收敛,极点全部在z = 0处(因果系统),因而只能用较高的阶数达到高的选择性。FIR数字滤波器的幅频特性精度较之于IIR数字滤波器低,但是线性相位,就是不同频率分量的信号经过fir滤波器后他们的时间差不变,这是很好的性质。FIR数字滤波器是有限的单位响应也有利于对数字信号的处理,便于编程,用于计算的时延也小,这对实时的信号处理很重要。 FIR滤波器因具有系统稳定,易实现相位控制,允许设计多通带(或多阻带)滤波器等优点收到人们的青睐。

IIR滤波器采用递归型结构,即结构上带有反馈环路。IIR滤波器运算结构通常由延时、乘以系数和相加等基本运算组成,可以组合成直接型、正准型、级联型、并联型四种结构形式,都具有反馈回路。同时,IIR数字滤波器在设计上可以借助成熟的模拟滤波器的成果,如巴特沃斯、契比雪夫和椭圆滤波器等,有现成的设计数据或图表可查,在设计一个IIR数字滤波器时,我们根据指标先写出模拟滤波器的公式,然后通过一定的变换,将模拟滤波器的公式转换成数字滤波器的公式。

滤波器的设计可以通过软件或设计专用的硬件两种方式来实现。随着MATLAB软件及信号处理工具箱的不断完善,MATLAB很快成为应用学科等领域不可或缺的基础软件。它可以快速有效地实现数字滤波器的设计、分析和仿真,极大地减轻了工作量,有利于滤波器设计的最优化。

2.1 IIR滤波器原理

冲激响应不变法是使数字滤波器在时域上模拟滤波器,但是它们的缺点是产生频率响应的混叠失真,这是由于从s平面到z平面是多值的映射关系所造成的。

双线性变换法是使数字滤波器的频率响应与模拟滤波器的频率响应相似的一种变换方法。为了克服多值映射这一缺点,我们首先把整个s平面压缩变换到某一中介的s1平面的一条横带里,再通过变换关系将此横带变换到整个z平面上去,这样就使得s平面与z平面是一一对应的关系,消除了多值变换性,也就消除了频谱混叠现象。

双线性法设计IIR数字滤波器的步骤:

1)将数字滤波器的频率指标{

v2-ec61fd4430db8a3415fe89ef2f6bba3b_b.jpg

}由Wk=(2/T)*tan(

v2-074f912b620a77ea2debcd0eea34a37a_b.jpg

),转换为模拟滤波器的频率指标{

v2-ec61fd4430db8a3415fe89ef2f6bba3b_b.jpg

}.

2)由模拟滤波器的指标设计H(s).

3)由H(s)转换为H(z)


v2-15f9d1ca327d0fb7f9090e56271edd26_b.jpg


2.2 FIR滤波器原理

FIR滤波器与IIR滤波器特点不同,IIR滤波器的相位是非线性的,若需线性相位则要采用全通网络进行相位校正。而有限长单位冲激响应(FIR)数字滤波器就可以做成具有严格的线性相位,同时又可以具有任意的幅度特性。

由于FIR系统的冲激响应就是其系统函数各次项的系数,所以设计FIR滤波器的方法之一可以从时域出发,截取有限长的一段冲激响应作为H(z)的系数,冲激响应长度N就是系统函数H(z)的阶数。只要N足够长,截取的方法合理,总能满足频域的要求。这种时域设计、频域检验的方法一般要反复几个回合,不像IIR DF设计靠解析公式一次计算成功。给出的理想滤波器频率响应是

v2-e7443a044ef8bee2168ad8427c210147_b.jpg

,它是w的周期函数,周期

v2-e8a5aa8367498fda4c5104836c6485dc_b.jpg

,因此可展开成傅氏级数

v2-811648329de7e51d5dfdefb02b1e63e5_b.jpg

由傅立叶反变换导出,即

v2-75c80922207e0fdcc65d24e8aa0e8429_b.jpg

,再将

v2-39a8304b56ca76d2e10ba25950b26e47_b.jpg

与窗函数

v2-fcd32a7001a9bed9640a23f3b72b736f_b.jpg

相乘就可以得到

v2-1578fda437a0e987fc17d709461608b0_b.jpg

v2-d4578a105e8b60a2c8c0380d5490dd12_b.jpg

v2-2245478dbe44e2b6006480d681f9a34f_b.jpg

v2-47356e37b7afd0a81e04af6deed1936f_b.jpg

的计算可采用傅氏变换的现成公式和程序,窗函数

v2-6eaf6aebb0cfd00d8844eeeda7edcaf5_b.jpg

也是现成的。但整个设计过程不能一次完成,因为窗口类型和大小的选择没有解析公式可一次算,整个设计可用计算机编程来做。

窗函数的傅式变换W(ejω)的主瓣决定了H(ejω)过渡带宽。W(ejω)的旁瓣大小和多少决定了H(ejω)在通带和阻带范围内波动幅度,常用的几种窗函数有:

矩形窗 w(n)=RN(n);

Hanning窗

v2-a88367734af7b5329f79a8eefce03c5a_b.jpg

Hamming窗

v2-ffc2631869d5007a764070240b36a5f8_b.jpg

Blackmen窗

v2-471121d3ef8302edc9f10a7ce475e618_b.jpg

Kaiser窗

v2-c4f08fa7da93a49061038f176107bf5d_b.jpg

式中Io(x)为零阶贝塞尔函数。

根据题目要求选择窗函数法设计FIR滤波器。

窗函数设计线性相位FIR滤波器步骤如下:

  1. 确定数字滤波器的性能要求,临界频率

v2-5a4f79696c653b5a2f205bd5679f1312_b.jpg

,滤波器单位脉冲响应长度N

    1. 根据性能要求,合理选择单位脉冲响应h(n)的奇偶对称性,从而确定理想频率响应

v2-8c9217064d2be92b5672014fde85ae3d_b.jpg

的幅频特性和相频特性

    1. 求理想单位脉冲响应

v2-b3c3981c36f6c77704748d09c2a5f7a8_b.jpg

,在实际计算中,可对

v2-8c9217064d2be92b5672014fde85ae3d_b.jpg

采样,并对其求IDFT的

v2-ade1a063e3dbca23c43a394f6e9f2c77_b.jpg

,用

v2-3eacafc075639c656800963152255633_b.jpg

代替

v2-b3c3981c36f6c77704748d09c2a5f7a8_b.jpg


    1. 选择适当的窗函数w(n),根据

v2-12dffe561b8032ffb730855cc69e0b35_b.jpg

求所需设计的FIR滤波器单位脉冲响应

v2-8c9217064d2be92b5672014fde85ae3d_b.jpg

,分析其幅频特性,若不满足要求,可适当改变窗函数形式或长度N,重复上述设计过程,以得到满意的结果


3.语音信号滤波

用自己设计的各滤波器分别对加噪的语音信号进行滤波,在Matlab中,FIR滤波器利用函数fftfilt对信号进行滤波,IIR滤波器利用函数filter对信号进行滤波。

函数fftfilt用的是重叠相加法实现线性卷积的计算。调用格式为:y=fftfilter(h,x,M)。其中,h是系统单位冲击响应向量;x是输入序列向量;y是系统的输出序列向量;M是有用户选择的输入序列的分段长度,缺省时,默认的输入向量的重长度M=512。

函数filter的调用格式:yn=filter(B,A.xn),它是按照直线型结构实现对xn的滤波。其中xn是输入信号向量,yn输出信号向量。

4.GUI设计系统界面

图形用户界面(graphical user interfaces ,GUI)则是由窗口、光标、按键、菜单、文字说明等对象(objects)构成的一个用户界面。用户通过一定的方法(如鼠标或键盘)选择、激活这些图形对象,使计算机产生某种动作或变化,比如实现计算、绘图等。在MATLAB中GUI是一中包含多种对象的图形窗口,并为GUI开发提供一个方便高效的集成开发环境GUIDE。GUIDE主要是一个界面设计工具集,MAYLAB将所有GUI支持度控件都集成在这个环境中,并提供界面外观、属性和行为响应方式的设置方法。GUIDE将设计好的GUI保存在一个FIG文件中,同时生成M文件框架。

FIG文件:包括GUI图形窗口及其所有后裔的完全描述,包括所有对象属性的属性值。它是一个二进制文件调用hsave课保存图形窗口时将生车该文件。M文件包括GUI设计、控件函数以及定义为子函数的用户控件回调函数,主要用于控制GUI展开时的各种特征。 GUI创建包括界面设计和控件编程两部分,主要步骤如下。第一步:通过设置GUIDE应用程序的选项来运行GUIDE;第二步:使用界面设计编辑器进行面设计;第三步:编写控件行为响应控制(即回调函数)代码。


四.设计步骤与结果

1、录音、保存与回放(以男音为例)

在Matlab软件平台下,利用wavread对语音进行采样,记住采样频率和采样点数。

代码如下:

t=5;

fs=11025;

xn=wavrecord(t*Fs,Fs,1);

plot(xn,'r')

wavwrite(xn,Fs,'d:\nanyin2')

wavplay(xn,fs)

2、语音信号的频谱

可以用函数fft对信号进行快速傅里叶变换,得到信号的频谱特性。代码如下:x=wavread('f:nanyin2.wav',5*fs);

N=length(x);

t=(1:N)/Fs;

f0=Fs/N;

n=1:N/2;

figure(1);

subplot(3,1,1);

plot(t, x);

title('男音信号时域波形');

xlabel('时间t');

ylabel('振幅A');

grid ;

y=fft(x,N);

subplot(3,1,2);

plot(k,abs(y(n)));

title('FFT变换后的频谱特性');

xlabel('频率/Hz');

ylabel('幅值/A');

grid;

subplot(3,1,3);

if y~=0

plot(k,20*log10(abs(y(n))));

end

xlabel('Hz');

ylabel('振幅/db');

title('FFT变换后声音的频谱特性');

grid ;

下图给出所录语音信号的时域波形和频谱。


v2-f4143707971288f5c871fb04e5ad7483_b.jpg


从图上分析可得到,男生的声音信号频率基本上在0—1000hz之间,主要集中在50到550低频部分。

3、设计数字信号滤波器,画出频率响应

根据分析所得的数字信号滤波器的频谱特点,按相关滤波器指标进行滤波器设计。

  1. 双线性变换法设计IIR数字信号滤波器

用到的函数及其调用方法如下:

[Bz, Az]=bilinear(B, A, FS) 其中,B, A分别是H(s)的分子多项式、分母多项式的系数向量,Bz, Az分别是H(z) 的分子多项式、分母多项式的系数向量,FS是抽样频率。

[N,Wn]=buttord(wp,ws,Ap,As,'s'); 其中,wp,ws分别是通带和阻带的截止频率,是归一化频率,实验的wp,ws是标量,Ap,As是通带和阻带的衰减,N是滤波器的阶数。

[z,p,k]=buttap(N) 其中,N是欲设计的低通原型滤波器的阶数,z,p,k是G(p)的零极点和增益。

[b,a]=lp2lp(b,a,Wn); 其中,b,a是G(p)的分子分母多项式的系数向量

给定滤波器参数fp=1000;fs=1800;Ap=1;As=50;,双线性设计IIR滤波器

实现滤波器设计的代码如下:

fp=1000;fs=1800;FS=11025;

wp=2*pi*fp/FS;ws=2*pi*fs/FS;

Ap=1;As=50;

wp1=2*FS*tan(wp/2);ws1=2*FS*tan(ws/2);

[N1,Wn]=buttord(wp1,ws1,Ap,As,'s');

[z,p,k]=buttap(N1);

[b,a]=zp2tf(z,p,k);

[bt,at]=lp2lp(b,a,Wn);

[bz,az]=bilinear(bt,at,FS);

freqz(bz,az,N1,FS);

运行程序,可以得到双线性变换法数字滤波器的频谱图:


v2-8a165c980a20505c48b172821bc1a6c8_b.jpg


当频率为为1000hz时,衰减很少,而当频率为1800hz时,其衰减为50dB,满足参数要求.

(2)窗函数法设计FIR数字信号滤波器

用到的函数及其调用方法如下:

fp=1000;fs=2000;FS=11025;

wp=2*pi*fp/FS;ws=2*pi*fs/FS;

deltaw=ws-wp;

M=ceil(6.6*pi/deltaw);

N=M+mod(M+1,2);

wdham=(hamming(N))';

hd=ideallp(wc,N);

h=hd.*wdham;

[db,mag,pha,grd,w]=myfreqz(h,[1]);

subplot(121),

plot(w/pi,mag);

legend('哈明窗');

subplot(122),

plot(w/pi,db);

legend('哈明窗');

运行程序,可以得到窗函数设计的法数字滤波器的频谱图:


v2-7e3c34744008f7cdcab4e9254f848ac4_b.jpg


观察第一幅图:当频率为1000hz(上图0.18)时,衰减为1dB,当频率为2000hz(上图0.36)时,衰减为50dB,满足参数要求。

观察第二幅图:当频率小于2000hz时,相位变化是呈线性的。

4、用设计的滤波器进行滤波

用设计的滤波器对带噪信号进行滤波,比较滤波前后语音信号的波形及频谱,并在一个窗口上显示出来。Matlab中,FIR滤波器利用函数ffifilt进行滤波,IIR滤波器利用函数filter进行滤波。


v2-a8007c9098cff0a6af6a5780cca9ec73_b.jpg



v2-4267f7f3abb43a428961d9272bf8b28e_b.jpg


由图像可以看出,所设计的滤波器基本可以满足实验要求,能够滤除信号中的高频成分,但是通过播放滤除后的声音信号可以看出,滤波后的有效信号的强度明显有所下降.

5、加噪音信号并设计合适滤波器滤除噪音信号

对原始语音信号加噪声利用MATLAB中的随机函数(rand或randn)产生噪声加入到语音信号中,模仿语音信号被污染,并对其频谱分析。Randn函数有两种基本调用格式:Randn(n)和Randn(m,n),前者产生n×n服从标准高斯分布的随机数矩阵,后者产生m×n的随机数矩阵。在这里,我们选用Randn(m,n)函数。语音信号添加噪声及其频谱分析的主要程序如下:

x=wavread('d:\nanyin2.wav',5*fs);

n= length (x) ; %求出语音信号的长度

Noise=0.01*randn(n,2); %随机函数产生噪声

Si=x+Noise; %语音信号加入噪声

sound(Si);

subplot(2,1,1);

plot(Si);title('加噪语音信号的时域波形');

S=fft(Si); %傅里叶变换

subplot(2,1,2);

plot(abs(S));

title('加噪语音信号的频域波形');


运行程序,结果如图所示。




v2-47edd4b8f1c33ed0af7f1b3ce0716ce9_b.jpg



从时域图上可以看到:噪声贯穿时域的整个过程。播放加噪语音时,无法辨别所录男音。

从频谱图可以看出,原始信号的主要频率集中在0Hz附近的低频部分,加上噪声信号后,在5000Hz两侧有个幅值很大的高频成分,即以上所加的噪声。

利用上面设计的FIR滤波器(fs改为1800)对加噪男音滤波如下:


v2-0b9e5b44c1f06b68d2fa423977eaae84_b.jpg



v2-e9f5b8ace1a403af4ad1245c9ed55f0e_b.jpg


由图可以看出噪声很好被滤除,男音很好的被还原。

6、男声、女声的比较

(1)用同样的方法录制一段女生的语音信号,做出时域波形和频谱图进行比较,如下图:


v2-9b1932c95c9ed37c8654ef70da974dd9_b.jpg


比较男声和女声之间的频谱差异:

通过比较可以看出,男生和女生的声音的频带范围还是由区别的,一般说来,男生的频率要低一些集中在50HZ到550HZ,女生的频率要高一些,集中在200HZ到700HZ,这与实际生活的经验是一致的。男音低沉,女音清脆。

7、GUI演示界面的设计


四.收获与体会

  1. 滤波器指标的选择:

在刚开始做题时,首先碰到的问题是滤波器的指标选择问题。通过上网查资料,我们得到人类的声音范围在20到20000赫兹,男生的声音范围是60到550HZ,女生在150到1200HZ,根据资料采用如下指标:

  1. 通带频率fp=1000hz
  2. 阻带频率fs=2000或1800hz
  3. 通带衰减Ap=1dB
  4. 阻带衰减As=50dB
  5. 采样频率FS=11025hz

2、所设计的IIR滤波器与FIR滤波器的频率特性的区别

从性能上来说,IIR滤波器传函数的极点可位于单位圆内的任何地方,因此能用较低的阶数获得高的选择性,所用的存储单元少,经济效率高。但它是以牺牲相位的线性为代价的。选择性越好,则相位非线性越严重。相反,FIR滤波器却可以得到严格的线性相位,然而由于FIR滤波器传输函数的极点固定在原点,所以只能用较高的阶数达到较高的选择性;对于同样的滤波器设计指标,FIR滤波器所要求的阶数可以比IIR滤波器高5到10倍,结果成本较高,信号延时较大。

从总体上说,IIR滤波器达到同样效果阶数少,延迟小,但有稳定性问题和非线性相位;FIR滤波器没有稳定性问题,是线性相位,但阶数多,延迟大。

3、IIR滤波器设计中遇到的问题及解决方法

问题:过渡带过窄,导致IIR滤波器非线性严重。

对于此问题的出现,是由于采样频率过大而过渡带过小引起的。刚开始设计IIR 滤波器时,采样频率为44100hz,过渡带为800hz,IIR滤波器的非线性非常严重,且阶数很高。为了解决这个问题,在录音的时候将采样频率降低,比如由原来的44100hz减少为11025hz.

4、问题:理想滤波器的频响是怎样的?

对于理想滤波器,在通带内纹波小,在过渡带内能够迅速衰减且呈线性。

实验的特色部分就是从噪音中滤出人的声音,探讨怎样能得到比较满意的效果

5、有用信号的加强

当一个信号经过滤波器后,有用信号的幅值会相对减少,因此要对信号进行加强。采取的方法是在画图形的时候在幅值前面乘以一个系数。

6、实验分工:

在此次实验中,主要负责程序的编程与调试,主要负责GUI界面的设计修改。

六、参考文献

1.宋宇飞,《数字信号处理实验与学习指导》,清华大学出版社2012.8

2.陈怀琛,《数字信号处理教程-matlab释义与实现》,电子工业出版社,2008.10

3 陈垚光,《精通MATLAB GUI设计》,电子工业出版社,2008.2

4.胡光书,《数字信号处理》,清华大学出版社,2011.10

  • 1
    点赞
  • 1
    评论
  • 11
    收藏
  • 打赏
    打赏
  • 扫一扫,分享海报

评论 1 您还未登录,请先 登录 后发表或查看评论
©️2022 CSDN 皮肤主题:深蓝海洋 设计师:CSDN官方博客 返回首页

打赏作者

sunny_chenxi

你的鼓励将是我创作的最大动力

¥2 ¥4 ¥6 ¥10 ¥20
输入1-500的整数
余额支付 (余额:-- )
扫码支付
扫码支付:¥2
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、C币套餐、付费专栏及课程。

余额充值