使用脉冲响应不变法设计切比雪夫I型IIR滤波器.doc

语音信号滤波去噪脉冲响应不变法设计的并联型切比雪夫I型滤波器学生姓名： 指导老师： 摘 要 本课程设计主要内容是使用脉冲响应不变法设计的并联型切比雪夫I型IIR滤波器，对一段含噪语音信号进行滤波去噪处理并根据滤波前后的波形和频谱分析滤波性能。本课程设计仿真平台为MATLAB7.0，开发工具是M语言编程。首先利用录音工具录制一段语音信号，并人为加入一单频噪声，然后对信号进行频谱分析以确定所加噪声频率，并设计滤波器进行滤波去噪处理，最后比较滤波前后的波形和频谱并进行分析。由分析结果可知，滤波器后的语音信号与原始信号基本一致，即设计的切比雪夫I型IIR滤波器能够去除信号中所加单频噪声，达到了设计目的。 关键词 滤波去噪；脉冲响应不变法；切比雪夫I型；IIR滤波器；MATLAB 1 引 言本课程设计主要解决在含噪情况下对语音信号的滤波去噪处理，处理时采用的是利用脉冲响应不变法设计的切比雪夫I型IIR滤波器。1.1 课程设计目的（1）熟悉使用MATLAB语言来解决一些简单的课程问题。（2）锻炼动手能力及独立思考并解决问题的能力（3）增加对书本知识的理解以及掌握，如：IIR滤波器的工作原理，切比雪夫滤波器的设计，脉冲响应法的应用等等；（4）对录制好的音频行加噪和过滤，比较前后语音变化；1.2 课程设计内容（1）使用录音软件；录制一段语音信号，要求为PCM编码格式 单声道8K 8位 7k/s的语音信号，绘制波形并观察频谱特点；（2）加上噪声并绘制加噪信号波形及观察其频谱特点；（3）设计滤波器，绘制出滤波器的频谱并检验滤波器是否合格；（4）最后使用设计好的滤波器对噪声信号进项过滤，绘制过滤后的波形及观察滤波前后频谱变化；（5）绘制滤波器的并联型结构图；2 设计原理录制一段语音信号，绘制波形并观察其频谱，给定相应技术指标，用脉冲响应不变法设计的一个满足指标的切比雪夫I 型IIR滤波器，对该语音信号进行滤波去噪处理，比较滤波前后的波形和频谱并进行分析。2.1 IIR滤波器IIR滤波器具有无限长脉冲响应，一般来说，所有的模拟滤波器都有无限长脉冲响应。因此，IIR滤波器设计的基本方法是利用复值映射将大家熟知的模拟滤波器变换为数字滤波器。这种IIR滤波器设计的基本方法存在途径之一2如下：实行滤波器变换SZ实行频带变换 S S设计模拟低通滤波器期望的IIR滤波器2.2 切比雪夫I型滤波器切比雪夫I型滤波器是电子滤波器的一种。它的特点是在通带具有等波纹响应。其振幅平方函数如2-1式所示： （2-1）式中，N为整数，称为滤波器的阶，是通带波纹因子（它与Rp有关）。其中N阶切比雪夫多项式如2-2式所示： （2-2）其中 切比雪夫滤波器存在两种可能的幅度平方响应的形状，当N分别为基数和偶数时如图2-1所示：图2-2 切比雪夫I型滤波器的幅度特性2.3 脉冲响应不变法脉冲响应不变法是一种将模拟滤波器转化为数字滤波器的基本方法。它通过对模拟滤波器的单位冲激响应h（t）等间隔抽样来获取对应数字滤波器的单位脉冲响应hk,即hk=h(t)|t=kT =h(kT)，其中T是等间隔抽样。不失一般性地，假定连续时间滤波器的系统函数为： （2-3）所对应的时域脉冲响应是： （2-4）对Tdhc(t)采样得到的离散时间滤波器的脉冲响应是： （2-5）离散时间滤波器的系统函数为： （2-6）3 设计步骤语音信号滤波去噪使用脉冲响应不变法设计的并联型切比雪夫IIR滤波器，设计流程如下图：图3-1 步骤图3.1 录制语音信号利用录音软件（如下图3-2）录制PCM编码格式 单声道8K 8位 7k/s的语音信号，语音内容不限；图3-2 录制语音信号软件界面3.2 绘制语音信号波形图使用MATLAB软件，将语音信号放到MATLAB的工作目录中；使用 x,fs.bits=wavread(ldd.wav)函数读入语音信号；源程序如下：x,fs,bits=wavread(ldd.wav); % 输入参数为文件的全路径和文件名，输出的第一个参数是每个样本的值，fs是生成该波形文件时的采样率，bits是波形文件每样本的编码位数。sound(x,fs,bits); % 按指定的采样率和每样本编码位数回放N=length(x); % 计算信号x的长度t=0:1/fs:(N-1)/fs; % 计算时间范围，样本数除以采样频率x=x; X=abs(fft(x); % 对原始信号和加噪信号进行fft变换,取幅度谱X=X(1:N/2); % 截取前半部分deltaf=fs/N; % 计算频谱的谱线间隔f=0:deltaf:fs/2-deltaf; % 计算频谱频率范围绘制的波形图如下：图3-3 原始语音信号的频谱图3.3 绘制加噪后信号波形图在原始语音信号上加入单频干扰信号，fn=3000Hz,其源程序如下：x,fs,bits=wavread(ldd.wav); % 输入参数为文件的全路径和文件名，输出的第一个参数是每个样本的值，fs是生成该波形文件时的采样率，bits是波形文件每样本的编码位数。sound(x,fs,bits); % 按指定的采样率和每样本编码位数回放N=length(x); % 计算信号x的长度fn=2000; % 单频噪声频率，此参数可改t=0:1/fs:(N-1)/fs; % 计算时间范围，样本数除以采样频率x=x; y=x+0.1*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：通过频谱图可以直观的看到，在3000Hz处有尖锐的单峰，与干扰频率一直；通过sound(y,fs,bits)回放该录音，可以明显的听到尖锐的单频噪鸣声；图3-4 加入噪声后的频谱对比图3.4 滤波器设计将数字滤波器的设计指标设为通带截止频率fn=3000,fp=fn-200=2800Hz,阻带频率fc=fn-50=2950HZ,通带波纹Rp=1，阻带波纹As=16dB；滤波器设计源程序如下：fs=8000;fn=3000;fp=fn-200;fc=fn-50; %定义通带和阻带截止频率Rp=1;As=16; % 定义通带波纹和阻带衰减 wp=(fp/fs)*2*pi; ws=(fc/fs)*2*pi; %计算对应的数字频率T=1; %定义采样间隔OmegaP=wp/T;OmegaS=ws/T; %截止频率线性变换Ripple=10(-Rp/20);Attn=10(-As/20); %转化为绝对形式cs,ds=afd_chb1(OmegaP,OmegaS,Rp,As); %计算滤波器阶数和截止频率b,a=imp_invr(cs,ds,T); % 脉冲响应不变法得到滤波器系统函数分子分母系数db,mag,pha,grd,w=freqz_m(b,a); %验证滤波器是否达到指定性能delta=1,zeros(1,99); ha=filter(b,a,delta); %计算脉冲响应运行程序得到的滤波器响应图如下：图3-5 滤波器响应图在命令行分别输入wp/pi,ws/pi验证滤波器是否满足设计标准： wp/pians = 0.7000 ws/pians = 0.73753.5 对噪声信号进行滤波，观察前后频谱变化使用设计好的滤波器对噪声信号进行过滤，其源程序如下：y_fil=filter(b,a,y); % IIR滤波器对信号y进行滤波处理，其中b,a为上面设计好的滤波器参数。Y_fil=abs(fft(y_fil); % 对fft变换,取幅度谱Y_fil=Y_fil(1:N/2); % 截取前半部分 sound(y_fil,fs,bits);程序运行结果如下图：图3-6 语音信号频谱分析图 由图3-6可以看出频率3000Hz被过滤掉，3000Hz处的单频尖峰噪声幅度明显减弱，使用sound(y_fil,fs,bits)函数播放过滤后的语音，可以听出过滤前后的语音变化，说明滤波器设计基本符合性能要求。3.6 绘制滤波器的并联型结构图设计完成滤波器后在窗体命令行输入C,B,A=dir2par(b,a)可以得出并联型结构的系数C B A，如下图所示：图3-7 并联型结构系数图根据并联型结构系数以及滤波器阶数可以画出滤波器的结构图如下：图3-8 滤波器并联型结构图4 异常处理（1）录制语音信号的时候不了解音频的编码方式的区别，所以最初录制的语音信号不能满足课设要求，后来使用3.1步骤中的录音软件解决了问题。（2）由于对MATLAB使用能力有限，在语音信号加噪时出了错，频谱图中没有出现单频噪声尖峰；然后通过网上搜索资料，并在同学的帮助下成功的显示了正确的频谱图；（3）在利用脉冲响应不变法设计切比雪夫滤波器时得不到理想的滤波器，把通带截止频率与阻带截止频率之间的差值设置的太小或太大，通过适当的调整参数，绘制出来的图像效果比较明显，基本符合设计指标。通过不断调整As，最后取As=16dB,这样通过MATLAB运算出来的滤波器的阻带波纹达到要求；（4）在设计切比雪夫I型IIR滤波器过程中，程序运行报错：Undefined command/afd_chb1，刚开始不知道怎么回事儿，后来老师告诉我此错误的原因在于没有将自编调用函数放进MATLAB的函数调用库中，将调用函数afd_chb1编写放入MATLAB的work文件夹后，问题就解决了；5 结束语课程设计的尘嚣渐渐远去，怀着些许不舍，怀着对当初课程设计开始时候的豪情万丈的决心的留恋，怀着通过这次课程设计积累的信心与斗志，我开始写下了这份结束语，为自己的足迹留下哪怕是微不足道但是对自己弥足珍贵的痕迹并期望与大家共勉。此次课程设计给我留下了太多记忆深刻的东西，学到了很多课本上学不到的知识。首先，我能够运用MATLAB，数字信号处理教程等学科所学的知识，并结合网络，图书馆搜集到的各种各样的材料，进行分析，比较，取舍，对课本知识有了一个深入的了解，这对我的学业生涯是大有裨益的。其次，这次课程设计使我有了一次通过理论联系实际，来解决实际问题的经历，它培养了自己分析问题，解决问题的能力，以及上网检索信息的能力；并且在和同学或者老师的交流中我体会到了团队的意义，以及探索未知学问的乐趣。其实学到的知识其实是次要的，重要的是我们探索知识的过程，这个过程便是一个人自主学习能力的体现，它将影响着我今后的发展。最后，这次的课程设计也激发了我对通行工程本专业的兴趣。曾经我对我们通信工程专业是一片迷茫，但是通过这次课程设计，通过亲自体验语音信号的加噪过滤过程，我似乎隐约了解了一些我们同形词专业的未来的发展方向。这种感悟，这种兴趣的来源只靠平时上课老师的说教是肯定不可能领悟到的。参考文献1 程佩青.数字信号处理教程M北京：清华大学出版社，20132 (美)维纳K英格尔 ,（美）约翰G普罗克斯（著）.刘树棠（译）.数字信号处理（MATLAB版）.第2版.西安：西安交通大学出版社，20083 张威. MATLAB基础与编程入门M.西安：西安电子科技大学出版社 2008附录1：语音信号滤波去噪设计源程序清单% 程序名称：finish.m% 程序功能：语音信号去噪使用脉冲响应不变法设计并联型切比雪夫I型滤波器% 程序作者：% 最后修改日期：x,fs,bits=wavread(ldd.wav); % 输入参数为文件的全路径和文件名，输出的第一个参数是每个样本的值，fs是生成该波形文件时的采样率，bits是波形文件每样本的编码位数。sound(x,fs,bits); % 按指定的采样率和每样本编码位数回放N=length(x); % 计算信号x的长度fn=2000; % 单频噪声频率，此参数可改t=0:1/fs:(N-1)/fs; % 计算时间范围，样本数除以采样频率x=x; y=x+0.1*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; % 计算频谱频率范围subplot(2,2,1);plot(t,x);title(原始信号);xlabel(时间/s);ylabel(幅度);axis(0,2.5,-1,1);grid onsubplot(2,2,2);plot(f,X);title(语音信号幅度谱);xlabel(频率,单位/Hz);ylabel(幅度谱);grid onsubplot(2,2,3);plot(t,y);title(加入一个单频信号的语音信号);xlabel(时间,单位/s);ylabel(幅度);axis(0,2.5,-1.1,1.1);grid onsubplot(2,2,4);plot(f,Y);title(加入一个单频信号的语音信号幅度谱);xlabel(频率,单位/Hz);ylabel(幅度谱);axis(0,4000,0,800);grid onfs=8000;fn=3000;fp=fn-200;fc=fn-50; %定义通带和阻带截止频率Rp=1;As=16; % 定义通带波纹和阻带衰减wp=fp/fs*2*pi;ws=fc/fs*2*pi; %计算对应的数字频率 T=1; %定义采样间隔OmegaP=wp/T; OmegaS=ws/T; %截止频率线性变换Ripple=10(-Rp/20);Attn=10(-As/20); %转化为绝对形式cs,ds=afd_chb1(OmegaP,OmegaS,Rp,As); %计算滤波器阶数和截止频率b,a=imp_invr(cs,ds,T);% 脉冲响应不变法得到数字滤波器系统函数分子分母系数db,mag,pha,grd,w=freqz_m(b,a); %验证滤波器是否达到指定性能delta=1,zeros(1,99);ha=filter(b,a,delta);%计算脉冲响应Subplot(221);plot(w/pi,db)%画出滤波器的幅度响应图xlabel(w/pi);ylabel(dB);title(滤波器幅度响应图)Subplot(222);plot(w/pi,mag)%画出滤波器幅度响应图xlabel(w/pi);ylabel(幅度mag);title(滤波器幅度响应)Subplot(223);plot(w/pi,pha)%画出滤波器的相位响应图xlabel(w/pi);ylabel(相位pha);title(滤波器相位响应)Subplot(224);plot(ha)%画出滤波器的脉冲响应图xlabel(n);ylabel(h(n);title(滤波器脉冲响应) y_fil=filter(b,a,y); % IIR