第15卷第1期2013年1月
大连民族学院学报JournalofDalianNationalitiesUniversity
Vol.15,No.1January2013
文章编号:1009-315X(2013)01-0072-04
巴特沃斯低通滤波器的实现方法研究
赵晓群,张
洁
(同济大学电子与信息工程学院,上海201804)
摘
要:从理论上分析巴特沃斯低通滤波器以及基于圆周卷积的设计方法,以截止频率为1kHz的巴特
沃斯低通滤波器为例,给出了一种基于C语言的通用巴特沃斯低通滤波器设计方法。从时间开销和处理与基于Matlab的设计方法进行比较。仿真结果表明:处理时间是基于Matlab的设计方法的时延等方面,1/200。
关键词:巴特沃斯滤波器;圆周卷积;C语言中图分类号:TN919
文献标志码:A
ResearchonImplementalmethodsofButterworthLow-PassFilter
ZHAOXiao-qun,ZHANGJie
(CollegeofElectronicsandInformationEngineering,TongjiUniversity,Shanghai201804,China)
Abstract:InthisarticleweanalyzedtheoreticallythedesignmethodsofButterworthlow-passfilterbasedoncircularconvolution.TakentheButterworthlow-passfilterwiththecutofffre-quency1kHzforanexample,wepresentadesignmethodofgeneralButterworthlow-passfil-terbasedonClanguage.Fromthetimeoverhead,processingdelayandsoon,wecomparethemethodbasedontheClanguagewiththemethodbasedonMatlab.Simulationresultsshowthattheprocessingtimeofthemethodgiveninthispaperis1/200ofthatofthemethodbasedonMatlab.
Keywords:Butterworthfilter;circularconvolution;Clanguage
巴特沃斯(Butterworth)滤波器是一种具有最大
平坦幅度响应的低通滤波器,它在通信领域里已有广泛应用。与贝塞尔(bessl)、契比雪夫(chebyshev)滤波器相比,巴特沃斯滤波器在线性相位、衰减斜率和加载特性3个方面具有特性均衡的优点,因此在实
[1]
际使用中,巴特沃斯滤波器已被列为首选。
Matlab语言是一种面向科学与工程计算的语它的测试程序手段丰富,扩展能力强,内涵丰言,
富。信号处理工具箱(SignalProcessingToolbox)提供了设计巴特沃斯滤波器的函数,本文利用这些函数,进行了巴特沃斯滤波器的程序设计,并将其作为函数文件保存,可方便地进行调用。虽然Matlab拥有强大的数值运算功能,但编程运算效率并不高,对运行内存要求巨大,通过本文可以证
Matlab对程序的处理时间比C明,在相同条件下,
增加了200倍。因此在工程应用中,常常将C与Matlab结合起来运用。
1巴特沃斯低通滤波器的特性
巴特沃斯低通滤波器的幅频特性模平方为
2
12
|Ha(jΩ)|=(,(1)
+(Ω/Ωc))N为滤波器的阶数,其中,Ωc为低通滤波器的截
1
|Ha(jΩ)|2=,止频率,当Ω=Ωc时,所以Ωc
2是滤波器的电压-3dB点或称半功率点。不同阶次N的巴特沃斯滤波器特性如图1,这一幅频特性具有如下特点:
收稿日期:2012-07-05;最后修回日期:2012-10-17
作者简介:赵晓群(1962-),男,黑龙江依安人,教授,博士,博士生导师,主要从事数字语音信号处理理论与应用研究。
第1期赵晓群,等:巴特沃斯低通滤波器的实现方法研究73
(1)最大平坦性:可以证明在Ω=0点,它的
前2N-1阶导数都等于零,这表示巴特沃斯滤波器在Ω=0附近一段范围内是非常平直的,它以
原点的最大平坦性来逼近理想低通滤波器,“最平响应”
即由此而得名。(2)通带、阻带下降的单调性:这种滤波器具有良好的相频特性。
(3)3dB的不变性:随着N的增加,频带边缘下降越陡峭,越接近理想特性。但不管N是多少,幅频特性都通过-3dB点。当Ω=Ωc时,特性以20NdB/dec速度下降[2]
。
图1
巴特沃斯滤波器幅频特性
2基于FFT的巴特沃斯低通滤波器设计
方法
在连续时间系统中,可以利用卷积的方法求系
统的零状态响应,这时,首先把激励信号分解为冲击函数序列,然后令每一冲激函数单独作用于系统求其冲激响应,最后把这些响应叠加即可得到系统对此激励信号的零状态响应。这个叠加的过程表
现为求卷积积分。在离散时间系统中,
可以采用大体相同的方法分析,由于离散信号本身就是一个不
连续的序列,因此,激励信号分解为脉冲序列的工作就很容易完成,对应每个样值的激励,系统得到对此样值的响应,每一响应也是一个离散时间序列,把这些序列叠加即得到零状态响应。
若长度为N1的序列x(n)与长度为N2的序
列h(n)做线卷积得到y(n)为有线长序列,其长度为N1+N2-1。共需要进行N1N2次乘法运算,
在N=N,需N2
12=N的情况下次乘法运算。
如果把求线性卷积改为求圆周卷积,并借助
FFT技术,有可能减少卷积所需的运算工作量。
如图2给出直接卷积与快速卷积两种方案原理
图。由图2(b)可见,快速卷积的过程中,共需要两次FFT,一次IFFT计算,相当于三次FFT的运算量。巴特沃斯低通滤波器的H(k)是知道的,数据已置于存储器中,实际只需要两个FFT运算量,如果假定N1=N2=N,经补零后点数为N1+N2-1≈2N,因而需要2×(Nlog22N)次复数乘法运算。此外,为完成X(k)与H(k)两序列相乘,还需要做2N次复乘。复乘运算全部次数为
2(Nlog22N)+2N=2N(1+log22N)。
(2)
显然,随着N值增大,式(2)的数字要比N2
显著减少
。
图2直接卷积与快速卷积原理方框图
当x(n)长度很长时,即N1N2,通常不允许等x(n)全部采集齐后再进行卷积,否则使输出相
对于输入有较长的延时,另外,若N1+N2-1太
大,h(n)要补上太多的零点,很不经济,且FFT的计算时间也要很长。为此,采用分段卷积的方法,即把x(n)分成长度与h(n)相仿的一段段,分别求出每段卷积的结果,然后用相应的方式把它们结合起来,便是总的输出。但是,N1可能很长,以至于趋于无穷大,以语音信号为例,如果不采用分段卷积的方法将迟迟不能给出结果,也无法找到
那样大的储备来满足N[2]
1的需要。但是即使采用分段卷积,也要根据语音的分帧,一帧帧输入语音,再对语音信号卷积,当前帧未输入完,不能进行卷积,因此也是有着相当大的时延,对于实时性要求较高的语音通信系统,这种方法没有本文给出的C语言实现的方法适用。
3巴特沃斯低通滤波器的MATLAB实现
Matlab的信号处理工具箱提供了有关巴特沃
斯滤波器的函数buttap,buttord,butter[3]
。
设计巴特沃斯滤波器的程序实现如下,Butter函数可在给定滤波器性能的情况下,选择巴特沃斯滤波器的阶数n和截止频率ωc,从而可利用
74
大连民族学院学报
第15卷
butter函数设计巴特沃斯滤波器的传递函数[n,
ωc]=buttord(ωp,ωs,Rp,Rs,'s')可得到满足性能的模拟巴特沃斯滤波器的最小阶数n及截止频率
ωc,其中ωp为通带的拐角频率,ωs为阻带的拐角频率,ωp和ωs的单位均为rad/s;Rs为通带区的
最大波动系数,
Rp为Rs阻带区的最小衰减系数,Rp和Rs的单位都为dB。[b,a]=butter(n,ωc,'s')可设计截止频率为ωc的n阶低通模拟巴特沃斯滤波器.利用buttord函数、
butter函数编制设计巴特沃斯低通滤波器的MATLAB函数文件butter-design.m,其清单如下:
Function[N,Wc,b,a]=butterdesign(Wp,Rp,Ws,As)
[N,Wc]=buttord(Wp,Ws,Rp,As,’s’);[b,a]=butter(N,Wc,’s’);[h,W]=freqs(b,a);
subplot(2,1,1);plot(W,abs(h));subplot(2,1,2);plot(W,angle(h));
在低速率语音编码中,基音周期的判别准确
性将直接影响解码端合成语音的质量,
MELP编码器的基音提取中,
首先要用截止频率为1kHz的6阶巴特沃斯低通滤波器对语音信号进行滤s(n)波,
消除语音帧中高频成分对基音周期估算的影响[4]
。在这里我们用Matlab实现用截止频率为1kHz的6阶巴特沃斯低通滤波器对语音信号s(n)来滤波。
采样的一组语音信号:e=[
5.4051,64.1490,80.6676,115.3568,48.4321,-35.8420,-90.2927,-74.3950,-97.3883,-72.8974]。调用函数[b,a]=butterπ(N,ωc),N为滤波器阶数6,ωc是归一化截止频率0.25。得到系数b=[
0.0011,0.0063,0.0158,0.0210,0.0158,0.0063,0.0011],a=[1.0000,-2.9785,4.1361,-3.2598,1.5173,-0.3911,0.0434]。根据系数就可以调用滤波器函数对信号实行滤波,y=filter(b,a,e),得到结果y
=[0.0057,0.1185,0.9043,3.9779,11.9714,26.8049,46.5178,63.2256,65.4751,45.0873]。y为滤波后的信号。
4巴特沃斯滤波器的C语言实现
利用重复计算由差分方程得出的递推公式来实现线性时不变离散时间系统,根据线性时不变系统的直接I型和直接II型或规范型结构,系统
的输入和输出满足如下差分方程[5]
:
y[n]-N
M
1
αky[n-k]=k=0
bkx[n-k]。
(3)
并具有如下系统函数:M
H(z)=
k=0
bkz
-k
M
(4)
1-k=1
akz
-k
M=N=6,将之前得到的系数a和b对应带入差分方程,在C语言下变成得到函数文件如下:main(){doublea[6]=
{2.9785,-4.1361,3.2598,-1.5173,0.3911,-0.0434};
doubleb[7]={0.0011,0.0063,0.0158,0.0210,
0.0158,0.0063,0.0011};
doublex[7]={0,0,0,0,0,0,0};doubley[7]={0,0,0,0,0,0,0};
doublee[10]={5.4051,64.1490,80.6676,115.3568,48.4321,-35.8420,-90.2927,-74.3950,-97.3883,-72.8974};doublep[10]={0,0,0,0,0,0,0,0,0,0};inti,j=0;
for(i=0;i<10;i++){x[6]=x[5];x[5]=x[4];x[4]=x[3];x[3]=x[2];x[2]=x[1];x[1]=x[0];x[0]=e[i];y[6]=y[5];y[5]=y[4];y[4]=y[3];y[3]=y[2];y[2]=y[1];y[1]=y[0];y[0]=b[0]*x[0]+b[1]*x[1]+b[2]*x[2]+b[3]*x[3]+b[4]*x[4]+b[5]*x[5]+b[
6]*x[6]+a[0]*y[1]+a[1]*y[2]+a[2]*y[3]+a[3]*y[4]+a[4]*y[5]+a[5]*y[6]
p[j++]=y[0];}for(j=0;j<=10;j++)printf(
调试运行得到结果依次为:0.005946,0.122325,0.918027,4.009936,12.023378,26.865400,46.562176,63.232224,65.438151,45.021945。
从结果可以看出,本文给出的基于C语言的设计方法与Matlab设计方法运行结果一致。
5运行时空开销的比较
可以用tic,
toc函数计算下在Matlab里面运行的时间为0.000303s。用clock函数对C函数计算得到的运行时间为0.000000333s。
将同样数据的10个信号扩展到30个,Matlab运行时间0.000394s,
C函数为0.00000133s。将同样的10个信号扩展到60个,
Matlab运行时
第1期赵晓群,等:巴特沃斯低通滤波器的实现方法研究75
间为0.000516s,
C函数为0.00000233s。将同样的信号扩展到150个,Matlab运行时间为0.001090s,
C函数运行时间为0.00000533s。C语言运行时间效率要比用Matlab高得多,对于数字语音信号这样的运算,一般信号数字运
算量是相当巨大的,当数字信号运算量越大,C语言进行计算的效率优势就越明显。
6查看变量使用空间
在Matlab中命令窗口,用whos和whosglobal
来查看所有局部和全局变量占用的内存大小[3]
,可以看到a和b是7为双精度数值各56字节,
e和y为输入和输出,是全部存在里面的,10个输入的话就有80字节,输出也是一样的。
在C函数中,除了输入输出共设立了4个数组,p数组只是为了显示用的,真正运行函数是不一定要用到的,
a数组有6位占用48字节,b,x,y有7位双精度值有56字节,但是C函数有最大的优势,就是它的输入输出占用的空间是可以变化的,这与Matlab必须先存在程序中不同,对于大量数值信号的计算来说,
这样是很费空间的,所以当计算数值数量越大时,
C函数优势越明显。7结论
本文给出了一种用C语言实现巴特沃斯低通
滤波器的方法,
并且通过与Matlab巴特沃斯滤波器一般方法以及另一种快速圆周卷积设计方法的
比较中,得出本文用的C语言实现的方法在时间开销以及运行时延等方面都比Matlab的方法要
好,
基于C语言的设计方法的处理时间是基于Matlab的设计方法的1/200。说明本文给出的算法在对于高实时性并且具有大输入信号量的通信系统中是具有很大优势的。不足之处由于大输入信号量的通信系统对于信号输入量要求很大,也只能分段传输处理,所以还是有一定的延时。本课题由水声通信与海洋信息技术教育部重点实验室(厦门大学)资助。参考文献:
1]李钟慎.基于MATLAB设计巴特沃斯低通滤波器
[J].信息技术,2003,27(3):49-50.
2]郑君里,应启行,杨为理.信号与系统[
M].北京:高等教育出版社,
2000.3]张雪英.数字语音处理及MATLAB仿真[M].北京:电
子工业出版社,
2010.4]韩笑蕾.(甚)低速率水声语音编码算法研究[D].上
海:同济大学,
2012:57-58.5]OPPENHEIMAV,SCHAFERRW,BUCKJR.Dis-crete-Timesignalprocessing[D].Cambridge:Massa-chusettsInstituteofTechnology,
1999:284-285.(责任编辑刘敏)
[[[[[
第15卷第1期2013年1月
大连民族学院学报JournalofDalianNationalitiesUniversity
Vol.15,No.1January2013
文章编号:1009-315X(2013)01-0072-04
巴特沃斯低通滤波器的实现方法研究
赵晓群,张
洁
(同济大学电子与信息工程学院,上海201804)
摘
要:从理论上分析巴特沃斯低通滤波器以及基于圆周卷积的设计方法,以截止频率为1kHz的巴特
沃斯低通滤波器为例,给出了一种基于C语言的通用巴特沃斯低通滤波器设计方法。从时间开销和处理与基于Matlab的设计方法进行比较。仿真结果表明:处理时间是基于Matlab的设计方法的时延等方面,1/200。
关键词:巴特沃斯滤波器;圆周卷积;C语言中图分类号:TN919
文献标志码:A
ResearchonImplementalmethodsofButterworthLow-PassFilter
ZHAOXiao-qun,ZHANGJie
(CollegeofElectronicsandInformationEngineering,TongjiUniversity,Shanghai201804,China)
Abstract:InthisarticleweanalyzedtheoreticallythedesignmethodsofButterworthlow-passfilterbasedoncircularconvolution.TakentheButterworthlow-passfilterwiththecutofffre-quency1kHzforanexample,wepresentadesignmethodofgeneralButterworthlow-passfil-terbasedonClanguage.Fromthetimeoverhead,processingdelayandsoon,wecomparethemethodbasedontheClanguagewiththemethodbasedonMatlab.Simulationresultsshowthattheprocessingtimeofthemethodgiveninthispaperis1/200ofthatofthemethodbasedonMatlab.
Keywords:Butterworthfilter;circularconvolution;Clanguage
巴特沃斯(Butterworth)滤波器是一种具有最大
平坦幅度响应的低通滤波器,它在通信领域里已有广泛应用。与贝塞尔(bessl)、契比雪夫(chebyshev)滤波器相比,巴特沃斯滤波器在线性相位、衰减斜率和加载特性3个方面具有特性均衡的优点,因此在实
[1]
际使用中,巴特沃斯滤波器已被列为首选。
Matlab语言是一种面向科学与工程计算的语它的测试程序手段丰富,扩展能力强,内涵丰言,
富。信号处理工具箱(SignalProcessingToolbox)提供了设计巴特沃斯滤波器的函数,本文利用这些函数,进行了巴特沃斯滤波器的程序设计,并将其作为函数文件保存,可方便地进行调用。虽然Matlab拥有强大的数值运算功能,但编程运算效率并不高,对运行内存要求巨大,通过本文可以证
Matlab对程序的处理时间比C明,在相同条件下,
增加了200倍。因此在工程应用中,常常将C与Matlab结合起来运用。
1巴特沃斯低通滤波器的特性
巴特沃斯低通滤波器的幅频特性模平方为
2
12
|Ha(jΩ)|=(,(1)
+(Ω/Ωc))N为滤波器的阶数,其中,Ωc为低通滤波器的截
1
|Ha(jΩ)|2=,止频率,当Ω=Ωc时,所以Ωc
2是滤波器的电压-3dB点或称半功率点。不同阶次N的巴特沃斯滤波器特性如图1,这一幅频特性具有如下特点:
收稿日期:2012-07-05;最后修回日期:2012-10-17
作者简介:赵晓群(1962-),男,黑龙江依安人,教授,博士,博士生导师,主要从事数字语音信号处理理论与应用研究。
第1期赵晓群,等:巴特沃斯低通滤波器的实现方法研究73
(1)最大平坦性:可以证明在Ω=0点,它的
前2N-1阶导数都等于零,这表示巴特沃斯滤波器在Ω=0附近一段范围内是非常平直的,它以
原点的最大平坦性来逼近理想低通滤波器,“最平响应”
即由此而得名。(2)通带、阻带下降的单调性:这种滤波器具有良好的相频特性。
(3)3dB的不变性:随着N的增加,频带边缘下降越陡峭,越接近理想特性。但不管N是多少,幅频特性都通过-3dB点。当Ω=Ωc时,特性以20NdB/dec速度下降[2]
。
图1
巴特沃斯滤波器幅频特性
2基于FFT的巴特沃斯低通滤波器设计
方法
在连续时间系统中,可以利用卷积的方法求系
统的零状态响应,这时,首先把激励信号分解为冲击函数序列,然后令每一冲激函数单独作用于系统求其冲激响应,最后把这些响应叠加即可得到系统对此激励信号的零状态响应。这个叠加的过程表
现为求卷积积分。在离散时间系统中,
可以采用大体相同的方法分析,由于离散信号本身就是一个不
连续的序列,因此,激励信号分解为脉冲序列的工作就很容易完成,对应每个样值的激励,系统得到对此样值的响应,每一响应也是一个离散时间序列,把这些序列叠加即得到零状态响应。
若长度为N1的序列x(n)与长度为N2的序
列h(n)做线卷积得到y(n)为有线长序列,其长度为N1+N2-1。共需要进行N1N2次乘法运算,
在N=N,需N2
12=N的情况下次乘法运算。
如果把求线性卷积改为求圆周卷积,并借助
FFT技术,有可能减少卷积所需的运算工作量。
如图2给出直接卷积与快速卷积两种方案原理
图。由图2(b)可见,快速卷积的过程中,共需要两次FFT,一次IFFT计算,相当于三次FFT的运算量。巴特沃斯低通滤波器的H(k)是知道的,数据已置于存储器中,实际只需要两个FFT运算量,如果假定N1=N2=N,经补零后点数为N1+N2-1≈2N,因而需要2×(Nlog22N)次复数乘法运算。此外,为完成X(k)与H(k)两序列相乘,还需要做2N次复乘。复乘运算全部次数为
2(Nlog22N)+2N=2N(1+log22N)。
(2)
显然,随着N值增大,式(2)的数字要比N2
显著减少
。
图2直接卷积与快速卷积原理方框图
当x(n)长度很长时,即N1N2,通常不允许等x(n)全部采集齐后再进行卷积,否则使输出相
对于输入有较长的延时,另外,若N1+N2-1太
大,h(n)要补上太多的零点,很不经济,且FFT的计算时间也要很长。为此,采用分段卷积的方法,即把x(n)分成长度与h(n)相仿的一段段,分别求出每段卷积的结果,然后用相应的方式把它们结合起来,便是总的输出。但是,N1可能很长,以至于趋于无穷大,以语音信号为例,如果不采用分段卷积的方法将迟迟不能给出结果,也无法找到
那样大的储备来满足N[2]
1的需要。但是即使采用分段卷积,也要根据语音的分帧,一帧帧输入语音,再对语音信号卷积,当前帧未输入完,不能进行卷积,因此也是有着相当大的时延,对于实时性要求较高的语音通信系统,这种方法没有本文给出的C语言实现的方法适用。
3巴特沃斯低通滤波器的MATLAB实现
Matlab的信号处理工具箱提供了有关巴特沃
斯滤波器的函数buttap,buttord,butter[3]
。
设计巴特沃斯滤波器的程序实现如下,Butter函数可在给定滤波器性能的情况下,选择巴特沃斯滤波器的阶数n和截止频率ωc,从而可利用
74
大连民族学院学报
第15卷
butter函数设计巴特沃斯滤波器的传递函数[n,
ωc]=buttord(ωp,ωs,Rp,Rs,'s')可得到满足性能的模拟巴特沃斯滤波器的最小阶数n及截止频率
ωc,其中ωp为通带的拐角频率,ωs为阻带的拐角频率,ωp和ωs的单位均为rad/s;Rs为通带区的
最大波动系数,
Rp为Rs阻带区的最小衰减系数,Rp和Rs的单位都为dB。[b,a]=butter(n,ωc,'s')可设计截止频率为ωc的n阶低通模拟巴特沃斯滤波器.利用buttord函数、
butter函数编制设计巴特沃斯低通滤波器的MATLAB函数文件butter-design.m,其清单如下:
Function[N,Wc,b,a]=butterdesign(Wp,Rp,Ws,As)
[N,Wc]=buttord(Wp,Ws,Rp,As,’s’);[b,a]=butter(N,Wc,’s’);[h,W]=freqs(b,a);
subplot(2,1,1);plot(W,abs(h));subplot(2,1,2);plot(W,angle(h));
在低速率语音编码中,基音周期的判别准确
性将直接影响解码端合成语音的质量,
MELP编码器的基音提取中,
首先要用截止频率为1kHz的6阶巴特沃斯低通滤波器对语音信号进行滤s(n)波,
消除语音帧中高频成分对基音周期估算的影响[4]
。在这里我们用Matlab实现用截止频率为1kHz的6阶巴特沃斯低通滤波器对语音信号s(n)来滤波。
采样的一组语音信号:e=[
5.4051,64.1490,80.6676,115.3568,48.4321,-35.8420,-90.2927,-74.3950,-97.3883,-72.8974]。调用函数[b,a]=butterπ(N,ωc),N为滤波器阶数6,ωc是归一化截止频率0.25。得到系数b=[
0.0011,0.0063,0.0158,0.0210,0.0158,0.0063,0.0011],a=[1.0000,-2.9785,4.1361,-3.2598,1.5173,-0.3911,0.0434]。根据系数就可以调用滤波器函数对信号实行滤波,y=filter(b,a,e),得到结果y
=[0.0057,0.1185,0.9043,3.9779,11.9714,26.8049,46.5178,63.2256,65.4751,45.0873]。y为滤波后的信号。
4巴特沃斯滤波器的C语言实现
利用重复计算由差分方程得出的递推公式来实现线性时不变离散时间系统,根据线性时不变系统的直接I型和直接II型或规范型结构,系统
的输入和输出满足如下差分方程[5]
:
y[n]-N
M
1
αky[n-k]=k=0
bkx[n-k]。
(3)
并具有如下系统函数:M
H(z)=
k=0
bkz
-k
M
(4)
1-k=1
akz
-k
M=N=6,将之前得到的系数a和b对应带入差分方程,在C语言下变成得到函数文件如下:main(){doublea[6]=
{2.9785,-4.1361,3.2598,-1.5173,0.3911,-0.0434};
doubleb[7]={0.0011,0.0063,0.0158,0.0210,
0.0158,0.0063,0.0011};
doublex[7]={0,0,0,0,0,0,0};doubley[7]={0,0,0,0,0,0,0};
doublee[10]={5.4051,64.1490,80.6676,115.3568,48.4321,-35.8420,-90.2927,-74.3950,-97.3883,-72.8974};doublep[10]={0,0,0,0,0,0,0,0,0,0};inti,j=0;
for(i=0;i<10;i++){x[6]=x[5];x[5]=x[4];x[4]=x[3];x[3]=x[2];x[2]=x[1];x[1]=x[0];x[0]=e[i];y[6]=y[5];y[5]=y[4];y[4]=y[3];y[3]=y[2];y[2]=y[1];y[1]=y[0];y[0]=b[0]*x[0]+b[1]*x[1]+b[2]*x[2]+b[3]*x[3]+b[4]*x[4]+b[5]*x[5]+b[
6]*x[6]+a[0]*y[1]+a[1]*y[2]+a[2]*y[3]+a[3]*y[4]+a[4]*y[5]+a[5]*y[6]
p[j++]=y[0];}for(j=0;j<=10;j++)printf(
调试运行得到结果依次为:0.005946,0.122325,0.918027,4.009936,12.023378,26.865400,46.562176,63.232224,65.438151,45.021945。
从结果可以看出,本文给出的基于C语言的设计方法与Matlab设计方法运行结果一致。
5运行时空开销的比较
可以用tic,
toc函数计算下在Matlab里面运行的时间为0.000303s。用clock函数对C函数计算得到的运行时间为0.000000333s。
将同样数据的10个信号扩展到30个,Matlab运行时间0.000394s,
C函数为0.00000133s。将同样的10个信号扩展到60个,
Matlab运行时
第1期赵晓群,等:巴特沃斯低通滤波器的实现方法研究75
间为0.000516s,
C函数为0.00000233s。将同样的信号扩展到150个,Matlab运行时间为0.001090s,
C函数运行时间为0.00000533s。C语言运行时间效率要比用Matlab高得多,对于数字语音信号这样的运算,一般信号数字运
算量是相当巨大的,当数字信号运算量越大,C语言进行计算的效率优势就越明显。
6查看变量使用空间
在Matlab中命令窗口,用whos和whosglobal
来查看所有局部和全局变量占用的内存大小[3]
,可以看到a和b是7为双精度数值各56字节,
e和y为输入和输出,是全部存在里面的,10个输入的话就有80字节,输出也是一样的。
在C函数中,除了输入输出共设立了4个数组,p数组只是为了显示用的,真正运行函数是不一定要用到的,
a数组有6位占用48字节,b,x,y有7位双精度值有56字节,但是C函数有最大的优势,就是它的输入输出占用的空间是可以变化的,这与Matlab必须先存在程序中不同,对于大量数值信号的计算来说,
这样是很费空间的,所以当计算数值数量越大时,
C函数优势越明显。7结论
本文给出了一种用C语言实现巴特沃斯低通
滤波器的方法,
并且通过与Matlab巴特沃斯滤波器一般方法以及另一种快速圆周卷积设计方法的
比较中,得出本文用的C语言实现的方法在时间开销以及运行时延等方面都比Matlab的方法要
好,
基于C语言的设计方法的处理时间是基于Matlab的设计方法的1/200。说明本文给出的算法在对于高实时性并且具有大输入信号量的通信系统中是具有很大优势的。不足之处由于大输入信号量的通信系统对于信号输入量要求很大,也只能分段传输处理,所以还是有一定的延时。本课题由水声通信与海洋信息技术教育部重点实验室(厦门大学)资助。参考文献:
1]李钟慎.基于MATLAB设计巴特沃斯低通滤波器
[J].信息技术,2003,27(3):49-50.
2]郑君里,应启行,杨为理.信号与系统[
M].北京:高等教育出版社,
2000.3]张雪英.数字语音处理及MATLAB仿真[M].北京:电
子工业出版社,
2010.4]韩笑蕾.(甚)低速率水声语音编码算法研究[D].上
海:同济大学,
2012:57-58.5]OPPENHEIMAV,SCHAFERRW,BUCKJR.Dis-crete-Timesignalprocessing[D].Cambridge:Massa-chusettsInstituteofTechnology,
1999:284-285.(责任编辑刘敏)
[[[[[