数组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设计流程图
图3.1总体设计流程图
图3.2最优等波纹法设计FIR滤波器流程图
本课程设计的总体设计流程如图3.1所示。
利用最优等波纹法设计FIR滤波器的流程图如图3.2所示。
3.2录制语音信号
在电脑上自己录制一段语音信号,且为单声道,再将此.
格式语音控制在3秒内,以减少设计中的误差。
然后在Matlab软件平台下,利用函数wavread对语音信号进行采样,记住采样频率和采样点数,采集完成后在信号中加入一个单频噪声,绘制原语音信号和加噪后的语音信号的时域和频域的波形图。
具体调用如下:
[x,fs,bits]=wavread('e:
\wangjun.wav');%输入参数为文件的全路径和文件名,输出的第一个参数是每个样本的值,fs是生成该波形文件时的采样率,bits是波形文件每样本的编码位数。
sound(x,fs,bits);%按指定的采样率和每样本编码位数回放
N=length(x);%计算信号x的长度
fn=2100;%单频噪声频率,此参数可改
t=0:
1/fs:
(N-1)/fs;%计算时间范围,样本数除以采样频率
x=x(:
1)';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.2所示:
图3-2原始语音信号与加噪后的语音信号时域、频域波形图
3.3滤波器设计
(1)设计数字滤波器
因为噪声的频率是1500Hz,所以设计各滤波器的性能指标如下所示:
,As=50dB,Rp=1dB
fp1=1400;fs1=1450;fs2=1550;fp2=1600;Rp=1;As=50;%带阻滤波器设计指标
wp1=fp1/Fs*2*pi;ws1=fs1/Fs*2*pi;ws2=fs2/Fs*2*pi;wp2=fp2/Fs*2*pi;%将Hz为单位的模拟频率换算为rad为单位的数字频率
[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
N=
288
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))
Asd=
43.7257
%校核最小阻带衰减并与As比较,然后增大(或降低)这个滤波器系数
%-------------
N=N+1;
h=firpm(N,f1,m,weights);
[db,mag,pha,grd,w]=freqz_m(h,[1]);%调用自编函数计算滤波器的频率特性
Asd=-max(db(ws1i:
1:
ws2i))
Asd=
42.3696
%校核最小阻带衰减并与As比较,然后增大(或降低)这个滤波器系数
%-------------中间省略了一些重复的步骤
N=N+1;
h=firpm(N,f1,m,weights);
[db,mag,pha,grd,w]=freqz_m(h,[1]);%调用自编函数计算滤波器的频率特性
Asd=-max(db(ws1i:
1:
ws2i))
Asd=
50.1366
%一直比较判断Asd和As的大小,当Asd第一次大于等于As时完成
%此时Asd>As,完成
N
N=
413
M=N+1
M=
414%滤波器长度M=N+1
(2)画图分析滤波器的各参数
n=0:
M;
subplot(2,2,1);stem(n,h);axis([0,M,-0.1,1])
xlabel('n');ylabel('h(n)');title('滤波器脉冲响应图');gridon
subplot(2,2,2);plot(w/pi,db);axis([0,1,-60,5]);
set(gca,'Xtick',[0,(wp1+ws1)/(2*pi),(wp2+ws2)/(2*pi),1],'Ytick',[-50,-1]);gridon
xlabel('w/pi');ylabel('dB');title('滤波器幅度响应图');
subplot(2,2,3);plot(w/pi,pha);axis([0,1,-4,4]);gridon
xlabel('w/pi');ylabel('相位');title('滤波器相位响应图');
subplot(2,2,4);plot(w/pi,mag);axis([0,1,0,1.1]);gridon
xlabel('w/pi');ylabel('幅度mag');title('滤波器幅度响应图');
图3.4最优等波纹的FIR滤波器
用最优等波纹法设计的FIR滤波器的各参数特性如图3.4所示。
观图可知,FIR滤波器符合设计要求。
3.4信号滤波处理
用滤波器对信号进行滤波
[y_fil,ny]=conv_m(y1,t,h,n);%用设计好的滤波器处理干扰后的音乐信号
y_fil=y_fil(1:
40000);%限定滤波后信号的长度
Y_fil=abs(fft(y_fil));%对原始信号和加噪信号进行fft变换,取幅度谱
Y_fil=Y_fil(1:
N/2);%截取前半部分
3.5结果分析
(1)在波形及频谱上比较滤波前后音乐信号
subplot(3,2,1);plot(t,y);axis([0,5,-1.3,1.3]);gridon
xlabel('时间(单位:
s)');ylabel('幅度');title('原始音乐信号');
subplot(3,2,2);plot(f,Y);gridon
xlabel('频率(单位:
Hz)');ylabel('幅度谱');title('音乐信号幅度谱图');
subplot(3,2,3);plot(t,y1);axis([0,5,-1.3,1.3]);gridon
xlabel('时间(单位:
s)');ylabel('幅度');title('加入单频干扰后的音乐信号');
subplot(3,2,4);plot(f,Y1);axis([0,4000,0,1500]);gridon
xlabel('频率(单位:
Hz)');ylabel('幅度谱');title('加入干扰后的音乐信号幅度谱图');
subplot(3,2,5);plot(t,y_fil);axis([0,5,-1.3,1.3]);gridon
xlabel('时间(单位:
s)');ylabel('幅度');title('滤波后的音乐信号');
subplot(3,2,6);plot(f,Y_fil);gridon
xlabel('频率(单位:
Hz)');ylabel('幅度谱');title('滤波后音乐信号幅度谱图');
图3.5滤波后信号幅度和频谱与原始信号和加噪后信号的比较
分析图3.5,通过FIR滤波器处理后的音乐信号和原始音乐信号在时域上基本相同,只是在频域1500Hz左右衰减较严重。
(2)在回放中比较滤波前后音乐信号
wavplay(y,Fs);%播放原始音乐信号
wavplay(y1,Fs);%播放加入单频噪声的音乐信号
wavplay(y_fil,Fs);%回放滤波后的音乐信号,与原始信号进行对比
通过在Matlab上回放原始音乐信号、加噪后的音乐信号和滤波后的音乐信号进行比较,滤波后的音乐信号在音质上和原始音乐基本相同,说明设计的FIR滤波器滤波效果较明显,设计成功。
4出现的问题及解决方法
1、加入的噪声如果不限定其振幅,对原始信号的影响较大。
为了便于观察和分析,限定其振幅为0.3的单频噪声(0.3*sin(fn*2*pi*t))。
2、在开始录制音乐信号并将其导入MATLAB中的时候出现过错误,原因是我录制的语音信号是双音频信号,不符合要求,在老师的指导下我将语音信号变成了单声道信号,再次导入的时候MATLAB不再报错。
3、在用最优等波纹法设计FIR滤波器的时候,因为本滤波器的设计指标相对较高,所以最后得到的滤波器的阶数会比较大,如果迭代运算时N的值每次只加1,会很耗时。
为了解决这个问题,在通过firpmord函数得到近似的滤波器的阶后,开始可以将N的值每次加大点的数值,到发现Asd的值和As很接近了的时候,再减小N每次加的数值,直到Asd刚大于或等于As。
4、在用设计好的FIR滤波器处理加噪信号时,不太会用fftfilt函数,所以就使用conv_m代替。
这两个函数的内部算法应该差不多,也达到了同样的效果。
5结束语
通过本次课程设计,我对数字信号处理有了更具体的认识,更深地了解了滤波器的设计一般设计步骤,尤其是对利用最优等波纹法设计FIR滤波器。
在本次课程设计中,我的题目是音乐信号的滤波去噪,相对平时作业题更贴近现实。
不像作业中设计滤波器最多是处理余弦信号等,处理音乐信号让我更形象地认识到滤波器在实际生活中的应用,也增强了学习兴趣。
还有就是在写课程设计报告时,你需要对原先课堂上学的知识有更透彻的认识,能对其原理和来龙去脉一清二楚。
所以写报告在一定程度上也是整理思路,加深对课程的理解。
课程设计是我们运用所学知识,动手实践的一个很好的机会。
它既可以帮助我们加深对所学知识的理解,又能提高我们运用知识,联系实际,动手实践的能力。
而且在设计过程中可能用到我们没学过的知识,需要我们去查阅资料获取相关信息,这又提高了我们查找信息和学习新知识的能力。
在实物的调试与检测过程中,又会遇到许多意想不到的问题,需要我们去分析原因和解决问题。
也体会到真正的去独立地完成一件事情是很困难的,同学以及老师的帮助和提醒是必须的。
通过这次课程设计,我拓宽了知识面,锻炼了实际操作能力,综合素质也得到了提高,进一步加深了了我们对专业的认识和激发了我们对专业的兴趣。
虽然课程设计结束了,但是我们的学习还没结束,对知识的进一步学习还需要继续,很开心成功地完成了这次设计。
6参考文献
[1][美]维纳·K·英格尔,约翰·G·普罗克斯.刘叔棠.数字信号处理(MATLAB版)DigitalSignalProcessingUsingMATLAB.第2版.西安:
西安交通大学出版社,2008年
[2]张志涌等.精通MATLAB6.5版.北京:
北京航空航天大学出版社,2003年
[3]邵玉斌.Matlab/Simulink通信系统建模与仿真实例分析.北京:
清华大学出版社,2008
[4]邓华.MATLAB通信仿真及应用实例详解.北京:
国防工业出版社,2003
附录1:
音乐信号滤波去噪——使用最优等波纹法设计的FIR滤波器
程序一:
%放出原音频
[x,fs,bits]=wavread('e:
\wangjun1.wav');
%输入参数为文件的全路径和文件名,
%输出的第一个参数是每个样本的值,
%fs是生成该波形文件时的采样率,
%bits是波形文件每样本的编码位数。
sound(x,fs,bits);%按指定的采样率和每样本编码位数回放
程序二:
%原音乐信号和加噪后的音乐信号的时域和频域的波形图绘制
N=length(x);%计算信号x的长度
fn=2100;%单频噪声频率,此参数可改
t=0:
1/fs:
(N-1)/fs;%计算时间范围,样本数除以采样频率
x=x(:
1)';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);
xlabel('时间(t)');ylabel('幅度');title('原始音乐信号');
axis([0,2.5,-1.5,1.5]);gridon;
subplot(2,2,2);plot(f,X);
xlabel('频率(f)');ylabel('幅度谱');title('原始音乐信号幅度谱');
axis([0,800,0,800]);gridon;
subplot(2,2,3);plot(t,y);
xlabel('时间(t)');ylabel('幅度');title('加干扰后的音乐信号');
axis([0,2.5,-1.5,1.5]);gridon;
subplot(2,2,4);plot(f,Y);
xlabel('频率(f)');ylabel('幅度谱');title('加干扰后的音乐信号幅度谱');
axis([0,800,0,800]);gridon;
程序三:
fp1=1400;fs1=1450;fs2=1550;fp2=1600;Rp=1;As=50;%带阻滤波器设计指标
wp1=fp1/Fs*2*pi;ws1=fs1/Fs*2*pi;ws2=fs2/Fs*2*pi;wp2=fp2/Fs*2*pi;%将Hz为单位的模拟频率换算为rad为单位的数字频率
[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
N=
288
h=firpm(N,f1,m,weights);%通过Parks-McClellan算法计算N阶FIR数字滤波器的频率响应
[db,mag,pha,grd,w]=freqz_m(h,[1]);%调用自编函数计