数组m包含有在f所给定频率上的期望幅度响应。
数组weights给出了再每个频带内的加权函数。
为了估计出滤波器的阶N,SP工具箱提供了函数firpmord。
这个函数也能估计在firpm函数中用到的其他参数。
这个基本语法是
[N,f0,m0,weights]=firpmord(f,m,delta);
这个函数计算窗的阶N,在f0中的归一化频带,m0中的振幅响应,以及在weights中的频带加权值。
向量f是归一化频带边缘频率向量,m是由f定义的频带上给定期望振幅值的向量。
F的长度小于两倍m的长度,也即f不包括0或1。
向量delta给定在每个频带内的容度。
估计出的参数现在能用于firpm函数。
3设计步骤
3.1设计流程图
本课程设计主要是对一段拉弦音乐信号,加入噪声后,用某种函数法设计出的FIR滤波器对加入噪声后的拉弦音乐信号进行滤波去噪处理,并且分析对比前后时域和频域波形的程序设计。
程序的设计流程图如下图所示:
图3-1总体设计流程图
图3-2最优等波纹法设计FIR滤波器流程图
本课程设计的总体设计流程如图3-1所示。
利用最优等波纹法设计FIR滤波器的流程图如图3-2所示。
3.2录制语音信号
音乐信号的采集
网上下载一段格式为.wav的音乐,注意调整格式为单声道,采样速率8kHz,8位码。
然后在Matlab软件平台下,利用函数wavread对语音信号进行采样,记住采样频率和采样点数。
通过wavread函数的使用,学生很快理解了采样频率、采样位数等概念。
采集完成后在信号中加入一个单频噪声,设计的任务即为从含噪信号中滤除单频噪声,还原原始信号。
[x,fs,bits]=wavread('lang.wav');%输入参数为文件的全路径和文件名,输出的第一个参数是每个样本的值,fs是生成该波形文件时的采样率,bits是波形文件每样本的编码位数。
sound(x,fs,bits);%按指定的采样率和每样本编码位数回放
N=length(x);%计算信号x的长度
fn=2100;%单频噪声频率,此参数可改
t=0:
1/fs:
(N-1)/fs;%计算时间范围,样本数除以采样频率
x=x';y=x+sin(fn*2*pi*t);
sound(y,fs,bits);%应该可以明显听出有尖锐的单频啸叫声
(1)原始信号和加噪信号的幅度和频谱分析
首先画出语音信号的时域波形;然后对语音号进行快速傅里叶变换,得到信号的频谱特性,从而加深对频谱特性的理解。
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-3原始信号和加入噪声后信号的幅度及其幅频特性图
原始信号和加噪信号的幅度和频谱如图3-3所示。
加入噪声后,原始信号幅度受到了一定的干扰。
加噪后的信号频谱图中也多了一个2100Hz的冲击频率。
3.3滤波器设计
设计数字滤波器和画出其频率响应
给出各滤波器的性能指标:
fp1=2000;fs1=2050;fs2=2150;fp2=2200;rp=1;as=50;
wp1=fp1/fs*2*pi;ws1=fs1/fs*2*pi;ws2=fs2/fs*2*pi;wp2=fp2/fs*2*pi;
[delta1,delta2]=db2delta(rp,as);%把以dB计的通带波纹和阻带衰减转换成相应的容度(或波纹)
f1=[wp1,ws1,ws2,wp2]/pi;m=[1,0,1];delta=[delta1,delta2,delta1];%设置firpmord的输入参数
[N,f1,m,weights]=firpmord(f1,m,delta);N%估猜窗的阶N
h=firpm(N,f1,m,weights);%通过Parks-McClellan算法计算N阶FIR数字滤波器的频率响应
[db,mag,pha,grd,w]=freqz_m(h,[1]);%调用自编函数计算滤波器的频率特性
delta_w=2*pi/1000;
ws1i=floor(ws1/delta_w)+1;wp1i=floor(wp1/delta_w)+1;
ws2i=floor(ws2/delta_w)+1;wp2i=floor(wp2/delta_w)+1;
asd=-max(db(ws1i:
1:
ws2i))
N=N+1;
h=firpm(N,f1,m,weights);
[db,mag,pha,grd,w]=freqz_m(h,[1]);%调用自编函数计算滤波器的频率特性
asd=-max(db(ws1i:
1:
ws2i))
N=N+1;
h=firpm(N,f1,m,weights);
[db,mag,pha,grd,w]=freqz_m(h,[1]);%调用自编函数计算滤波器的频率特性
asd=-max(db(ws1i:
1:
ws2i))
N
M=N+1
N=
288
asd=
42.1972
Warning:
OddordersymmetricFIRfiltersmusthaveagainofzero
attheNyquistfrequency.Theorderisbeingincreasedbyone.
Alternatively,youcanpassatrailing'h'argument,
asinfirpm(N,F,A,W,'h'),todesignatype4linearphasefilter.
>Insignal\private\firpminitat117
Infirpmat132
asd=
42.4740
asd=
42.4740
N=
290
M=
291
图3-4最优等波纹的FIR滤波器
用最优等波纹法设计的FIR滤波器的各参数特性如图3-4所示。
观图可知,FIR滤波器符合设计要求。
3.4信号滤波处理
在将滤波器设计好后,我们用自己设计的带阻滤波器对采集的语音信号进行滤波。
在Matlab中,FIR滤波器利用函数fftfilt对信号进行滤波。
我对信号进行滤波处理对应的程序如下:
[y_fil,ny]=conv_m(y,t,h,n);%用设计好的滤波器处理干扰后的音乐信号
y_fil=y_fil(1:
4000);%限定滤波后信号的长度
Y_fil=abs(fft(y_fil));%对原始信号和加噪信号进行fft变换,取幅度谱Y_fil=Y_fil(1:
N/2);%截取前半部分
在将加噪信号滤波之后,我们将滤波前后语音信号的波形及频谱图相互比较。
在同一张大图里分别绘制原始信号x,加噪信号y,滤波去噪信号y_fil的时域波形和频谱,以便比较和分析。
经过这段程序画出来的三个信号的时域波形和频谱图如下图3-5所示。
图3-5滤波后信号幅度和频谱与原始信号和加噪后信号的比较
分析图3-5,通过FIR滤波器处理后的原信号与滤波去噪信号的时域图基本相似;原信号与滤波去噪信号的频谱图波形也大致相似。
通过观察可以看到,加噪信号的时域图中大部分都被加入的噪声给遮盖了,加噪信号的频谱图中,我们可以很明显地看到与原信号频谱图相比,它在2100Hz左右处有一个尖脉冲,而滤波去噪信号的频谱图中
尖脉冲已经消失,波形大致与原图相似。
在将三个信号的时域波形和频谱图比较之后,我们还要通过回放去滤波去噪语音信号,来跟原信号相比,以检验滤波器的效果。
在Matlab中,函数sound可以对声音进行回放。
其调用格式为:
sound(x,fs,bits)。
我用sound(y_fil,fs,bits)语句回放该滤波去噪信号,便可以感觉到滤波后的语音信号与原信号差不多,但仍有一点点变化。
我们还需要将加噪后的信号和滤波之后的信号的频谱图进行分析比较。
图3-6加干扰与滤波后的信号幅度谱
3.5结果分析
在回放中比较滤波前后音乐信号
>>wavplay(y,Fs);%播放原始音乐信号
>>wavplay(y1,Fs);%播放加入单频噪声的音乐信号
>>wavplay(y_fil,Fs);%回放滤波后的音乐信号,与原始信号进行对比
通过在Matlab上回放原始音乐信号、加噪后的音乐信号和滤波后的音乐信号进行比较,可以感觉滤波前后的声音有变化,声音变得不会有刺耳的声音,而且比以前更加地平滑,滤波后的音乐信号在音质上和原始音乐基本相同,说明设计的FIR滤波器滤波效果较明显,设计成功。
4出现的问题及解决方法
1)在录音时,没有将录音的属性改为wav的格式,当在Matlab软件平台上调用时,出现无法调用的提示,最后通过反复思考终于找到问题的所在,在录音机的窗口上,点击“文件”,然后选项里点击“属性”中的“立即转换”可以将拉弦音乐的格式改为.wav的格式。
2)在用最优等波纹法设计FIR滤波器的时候,因为本滤波器的设计指标相对较高,所以最后得到的滤波器的阶数会比较大,如果迭代运算时N的值每次只加1,会很耗时。
为了解决这个问题,在通过firpmord函数得到近似的滤波器的阶后,开始可以将N的值每次加大点的数值,到发现Asd的值和As很接近了的时候,再减小N每次加的数值,直到Asd刚大于或等于As。
3)在用设计好的FIR滤波器处理加噪信号时,不太会用fftfilt函数,所以就使用conv_m代替。
这两个函数的内部算法应该差不多,也达到了同样的效果。
4)绘制出滤波后的波形,发现FIR滤波器没有滤掉单频噪声。
通过自己的仔细检查,是单频噪声的频率改动后,FIR滤波器的频率没有改动。
所以单频噪声的频率也应该自己先定义,FIR滤波器的截止频率应该以单频噪声的频率为中心,这样重新运行后,结果正确。
5结束语
通过本次课程设计,我对数字信号处理有了更具体的认识,更深地了解了滤波器的设计一般设计步骤,尤其是对利用最优等波纹法设计FIR滤波器。
这次课程设计比较的简单,老师给我们上了几节有关怎么样编写程序的课以及给了我们一些有关课程设计的模板,所以编写出来的程序就比较顺利地运行成功得出结果。
在本次课程设计中,我的题目是音乐信号的滤波去噪,相对平时作业题更贴近现实。
不像作业中设计滤波器最多是处理余弦信号等,处理音乐信号让我更形象地认识到滤波器在实际生活中的应用,也增强了学习兴趣。
还有就是在写课程设计报告时,你需要对原先课堂上学的知识有更透彻的认识,能对其原理和来龙去脉一清二楚。
所以写报告在一定程度上也是整理思路,加深对课程的理解,从这次课程设计中,我体会到什么事情都要自己真真正正用心的去做,才会使自己更加成长,没有学习就不可能有实践的机会,没有实践的机会就不会有所突破。
经过几天的设计和实践,已经比较顺利地完成了本次课程设计,在此感谢胡老师的指导和一些同学的帮助!
6参考文献
[1][美]维纳·K·英格尔,约翰·G·普罗克斯.刘叔棠.数字信号处理(MATLAB版)DigitalSignalProcessingUsingMATLAB.第2版.西安:
西安交通大学出版社,2008年
[2]张志涌等.精通MATLAB6.5版.北京:
北京航空航天大学出版社,2003年
[3]邓华.MATLAB通信仿真及应用实例详解.北京:
国防工业出版社,2003年
[4]程佩青.数字信号处理教程[m].北京:
清华大学出版社,2002年
[5]邵玉斌.Matlab/Simulink通信系统建模与仿真实例分析.北京:
清华大学出版社,2008年
[6]陈其宗.数字信号处理频谱计算与滤波器设计.北京:
电子工业出版社,2002
附录1:
用最优等波纹法设计FIR滤波器程序清单
%程序名称:
FIRLANG
%程序功能:
采用最优等波纹法设计FIR滤波器的方法,给弹拨音乐出去噪声。
%程序作者:
郎清凯
%最后修改日期:
2012-3-13
>>[x,fs,bits]=wavread('lang.wav');%输入参数为文件的全路径和文件名,输出的第一个参数是每个样本的值,fs是生成该波形文件时的采样率,bits是波形文件每样本的编码位数。
sound(x,fs,bits);%按指定的采样率和每样本编码位数回放
N=length(x);%计算信号x的长度
fn=2100;%单频噪声频率,此参数可改
t=0:
1/fs:
(N-1)/fs;%计算时间范围,样本数除以采样频率
x=x';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;%计算频谱频率范围
>>subplot(3,2,1);plot(t,x);axis([0,3,-2.3,2.3]);gridon;xlabel('时间(单位:
s)');ylabel('幅度');title('原始语音信号');
subplot(3,2,2);plot(f,X);gridon;xlabel('频率(单位:
Hz)');ylabel('幅度谱');title('语音信号幅度谱图');
subplot(3,2,3);plot(t,y);axis([0,3,-2.3,2.3]);gridon;xlabel('时间(单位:
s)');ylabel('幅度');title('加入单频干扰后的语音信号');
subplot(3,2,4);plot(f,Y);axis([0,4000,0,1500]);gridon;xlabel('频率(单位:
Hz)');ylabel('幅度谱');title('加入干扰后的语音信号幅度谱图');
>>fp1=2000;fs1=2050;fs2=2150;fp2=2200;rp=1;as=50;
wp1=fp1/fs*2*pi;ws1=fs1/fs*2*pi;ws2=fs2/fs*2*pi;wp2=fp2/fs*2*pi;
[delta1,delta2]=db2delta(rp,as);%把以dB计的通带波纹和阻带衰减转换成相应的容度(或波纹)
f1=[wp1,ws1,ws2,wp2]/pi;m=[1,0,1];delta=[delta1,delta2,delta1